Computer-implemented systems and methods for controlling sand production in a geomechanical reservoir system

ABSTRACT

Systems and methods are provided for use in predicting sand production in a geomechanical reservoir system. Computation of the sand production predictions can include solving a system of partial differential equations that model the geomechanical reservoir system. Systems and methods also are provided for use in operating a geomechanical reservoir system based on the sand production prediction for controlling sand production in the geomechanical reservoir system.

CROSS-REFERENCE TO RELATED APPLICATION

This application is a continuation of U.S. application Ser. No.12/561,886, filed Sep. 17, 2009, the entire contents of which areincorporated herein by reference.

1. TECHNICAL FIELD

This document relates to computer-implemented systems and methods foruse in predicting sand production in a geomechanical reservoir system.This document also relates to systems and methods for use in controllingsand production in a geomechanical reservoir system based on theprediction.

2. BACKGROUND

Production in a reservoir system is generally the phase that occursafter development of the reservoir, during which reservoir fluids, suchas hydrocarbons (oil or gas), are drained from the reservoir. Sanding isan occurrence in which formation solid particles are produced withreservoir fluids. The generic term used to describe small particles ofthe formation (the rock around a wellbore) which may be produced withthe reservoir fluid is “sand.” The term “fines” has been used in someliterature. Reservoir formation material generally comprises a rock typehaving sufficient porosity and permeability to store and transmitreservoir fluids, such as oil or water. Since sedimentary rocks areporous and form under temperature conditions at which fossil remnants(from which hydrocarbons are derived) can be preserved, they are themost common type of reservoir rocks (rather than igneous or metamorphicrocks). Examples of sedimentary rocks include, but are not limited to,conglomerates, sandstones, siltstones, shales, limestone, dolostone,halite (salt), salts, gypsum, and calcium sulfate anhydrite. Sedimentaryrocks can include a wide variety of minerals, including but not limitedto, quartz, feldspar, calcite, dolomite, and clay group minerals.

Sand production is the migration of the formation sand caused by theflow of reservoir fluids (such as oil) during production. Sandproduction can result from shear failure or tensile failure within thereservoir formation material. Shear failure can occur when the boreholepressure is significantly reduced (such as later in the life of a well),which increases stress near the wellbore, leading to formation failure.Tensile failure can occur when the porosity and permeability of thereservoir formation material near the wellbore are significantly damagedor when flow rates are extremely high. Under either tensile failurecondition, the flowing fluid can exert significant drag forces onindividual grains in the formation, which, if excessive, can cause thecementation between individual grains to fail, resulting in tensilefailure and sand production.

In many instances, sand production can be undesirable since it canrestrict productivity, erode completion components, impede wellboreaccess, interfere with the operation of downhole equipment, and presentsignificant disposal difficulties. Sand production increases thediameter of the perforation cavity or the borehole, reducing the supportaround the casing. As a result, perforations collapse, the cavitybecomes larger, and eventually production from the wellbore could cease.If sand production is severe, remedial action may be required to controlor prevent sand production altogether, such a gravel packing or sandconsolidation. In extreme cases, massive sanding can occur, in which thesand production increases uncontrollably, eventually completely erodingthe reservoir material forming the foundation of the well.

Conventionally, measures have been put into place early in a reservoirdevelopment project to prevent sand production altogether, to the extentpossible. For example, unconsolidated formations with multiple producingzones may be completed cased hole. The cased hole completion procedureinvolves constructing a barrier system around the borehole of thereservoir to prevent or retard the onset or extent of sanding, includingcementing steel pipes into the wellbore and using meshes, such as anexpandable sand screen technology. Such measures, however, may have anegative impact on well productivity. That is, productivity performanceof the cased hole well can be much lower than that of an open hole well.Also, such measures increase the cost of a reservoir developmentproject, since they require resources, time, and labor forimplementation.

The allowance of some amount of sand production can aid in increasingwell productivity. Therefore, a method for use in predicting sandproduction from a geomechanical reservoir system would be useful, inthat it would allow determination of if and when sand production wouldoccur, and to what extent, early on a well development project. Usingthese predictions, a reservoir can be constructed and operated such thata limited amount of sand is produced, a cavity generated in thereservoir is substantially sustained, and well productivity isincreased.

3. SUMMARY

As disclosed herein, computer-implemented systems and methods areprovided for use in predicting sand production from a geomechanicalreservoir system. The methods and systems comprise receiving dataindicative of gradual sanding and massive sanding of material within thegeomechanical reservoir system; computing, on a computer system, a valueof a critical plastic strain based on a fit of a geomechanical modelcomprising a hardening model to the received massive sanding data;computing, on a computer system, a value of at least one parameter of aproduction function based on a fit of a geomechanical model comprising aproduction function to the received gradual sanding data and using thevalue of the critical plastic strain and at least one value of one ormore hardening parameters; wherein the at least one value of the one ormore hardening parameters is computed based on a fit of the hardeningmodel to data indicative of plastic deformation of material within thegeomechanical reservoir system; wherein the hardening model modelsplastic deformation behavior of the material within the geomechanicalreservoir system; and wherein the production function predicts sandproduction from the geomechanical reservoir system. At least one valueof the one or more hardening parameters, the value of the criticalplastic strain, or the value of at least one parameter of the productionfunction may be output to a user interface device, a monitor, a computerreadable storage medium, a local computer, or a computer that is part ofa network.

Computer-implemented systems and methods also are provided for use inpredicting sand production from a geomechanical reservoir system,comprising: receiving data indicative of plastic deformation, gradualsanding and massive sanding of material within the geomechanicalreservoir system; computing, on a computer system, a value of one ormore hardening parameters based on a fit of a hardening model to thereceived plastic deformation data; wherein the hardening model modelsplastic deformation behavior of the material within the geomechanicalreservoir system; computing, on a computer system, a value of a criticalplastic strain based on a fit of a geomechanical model comprising thehardening model to the received massive sanding data; computing, on acomputer system, a value of at least one parameter of a productionfunction based on a fit of a geomechanical model comprising theproduction function to the received gradual sanding data and using avalue of at least one value of the hardening parameters and the value ofthe critical plastic strain; wherein the production function predictssand production from the geomechanical reservoir system. At least onevalue of the one or more hardening parameters, the value of the criticalplastic strain, or the value of at least one parameter of the productionfunction may be output to a user interface device, a monitor, a computerreadable storage medium, a local computer, or a computer that is part ofa network.

In one aspect of the foregoing methods and systems, the productionfunction is a function ƒ(x), where ƒ(x)=0 when x=0, and ƒ(x)=1 when x=1;and where x is a function of the critical plastic strain. In anotheraspect of the foregoing methods and systems, the value of at least oneparameter of the production function is a value of an exponent of theproduction function. The sand production function may be given by theexpression:

${f(x)} = \left( \frac{ɛ^{p}}{ɛ_{\lim}^{p}} \right)^{m}$

where x=ε^(p)/ε_(lim) ^(p), where ε^(p) is a plastic strain invariant,wherein ε_(lim) ^(p) is the critical plastic strain; and where m thevalue of the exponent.

The hardening model may be a modified Bigoni-Piccolroaz model. Theplastic deformation data may be obtained from a triaxial test. Thegradual sanding and massive sanding may be obtained from a hollowcylinder test. The fit of the hardening model to the received plasticdeformation data may be obtained by a regression.

In an aspect of the foregoing methods and systems, the fit of thegeomechanical model comprising the hardening model to the receivedmassive sanding data is obtained by solving a system of partialdifferential equations that model the geomechanical reservoir system;where the system of partial differential equations comprises a reservoirflow model and the geomechanical model comprising the hardening model;where the system of partial differential equations is coupled through afully-expanded Jacobian; and where the solving of the system of partialdifferential equations includes solving simultaneously in a single timestep the fully-expanded Jacobian based upon the received massive sandingdata.

In another aspect of the foregoing methods and systems, the fit of thegeomechanical model comprising the production function to the receivedgradual sanding data is obtained by solving a system of partialdifferential equations that model the geomechanical reservoir system;where the system of partial differential equations comprises a reservoirflow model and the geomechanical model comprising the hardening model;where the system of partial differential equations is coupled through afully-expanded Jacobian; and where the solving of the system of partialdifferential equations includes solving simultaneously in a single timestep the fully-expanded Jacobian based upon the received gradual sandingdata.

Computer-implemented systems and methods also are provided for use inpredicting sand production from a geomechanical reservoir system,comprising: receiving data indicative of physical or chemical propertiesassociated with material within the geomechanical reservoir system;generating sand production predictions by solving a system of partialdifferential equations that model the geomechanical reservoir system;wherein the system of partial differential equations comprises areservoir flow model and a geomechanical model of the geomechanicalreservoir system; wherein the geomechanical model comprises a hardeningmodel; wherein a sand production criterion is applied to thegeomechanical model; wherein the system of partial differential arecoupled through a fully-expanded Jacobian; wherein the solving of thesystem of partial differential equations includes solving simultaneouslyin a single time step the fully-expanded Jacobian based upon thereceived physical properties data; and wherein the generating isimplemented on a computer system. The generated sand productionpredictions may be output to a user interface device, a monitor, acomputer readable storage medium, a local computer, or a computer thatis part of a network. The sand production criterion may be a criticalvalue of a total strain, a critical value of a plastic strain invariant,or a maximum effective stress.

In addition, computer-implemented systems and methods also are providedfor use in predicting sand production in a geomechanical reservoirsystem, comprising: receiving data indicative of physical or chemicalproperties associated with material within the geomechanical reservoirsystem; defining a grid comprising a plurality of grid cells; generatingsand production predictions by solving a system of partial differentialequations that model the geomechanical reservoir system; wherein thesystem of partial differential equations comprises a reservoir flowmodel and a geomechanical model of the geomechanical reservoir system;wherein the geomechanical model comprises a hardening model; wherein asand production criterion is applied to the geomechanical model; whereinthe system of partial differential equations are coupled through afully-expanded Jacobian; wherein the solving of the system of partialdifferential equations includes solving simultaneously in a single timestep the fully-expanded Jacobian based upon the received physicalproperties data; wherein the reservoir model and the geomechanical modelare computed on the grid; and wherein the generating is implemented on acomputer system. The generated sand production predictions may be outputto a user interface device, a monitor, a computer readable storagemedium, a local computer, or a computer that is part of a network. Thesand production criterion may be a critical value of a total strain, acritical value of a plastic strain invariant, or a maximum effectivestress.

An aspect of the present disclosure provides a computer system forperforming the steps of any of the methods and systems disclosed herein,including the foregoing systems and methods. The computer systemcomprises one or more processor units; and one or more memory unitsconnected to the one or more processor units, the one or more memoryunits containing one or more modules which comprise one or more programswhich cause the one or more processor units to execute steps comprisingperforming the steps of any of the systems and methods disclosed herein,including the foregoing systems and methods. In the foregoingembodiments, the one or more memory units may contain one or moremodules which comprise one or more programs which cause the one or moreprocessor units to execute steps comprising outputting to a display, auser interface device, a tangible computer readable data storageproduct, or a tangible random access memory, a result of the systems andmethods. For example, as is applicable to the method being executed, theresult of the system or method which is output can be at least one valueof the one or more hardening parameters, the value of the criticalplastic strain, the value of at least one parameter of the productionfunction, or the generated sand production predictions.

Another aspect of the present disclosure provides a computer-readablemedium storing a computer program executable by a computer forperforming the steps of any of the systems and methods disclosed herein,including the foregoing systems and methods. A computer program productis provided for use in conjunction with a computer having one or morememory units and one or more processor units, the computer programproduct comprising a computer readable storage medium having a computerprogram mechanism encoded thereon, wherein the computer programmechanism can be loaded into the one or more memory units of thecomputer and cause the one or more processor units of the computer toexecute steps comprising performing the steps of any of the systems andmethods disclosed herein, including the foregoing systems and methods.In the foregoing embodiments, the computer program mechanism may beloaded into the one or more memory units of said computer and cause theone or more processor units of the computer to execute steps comprisingoutputting to a display, a user interface device, a tangible computerreadable data storage product, or a tangible random access memory, aresult of the system or method. For example, as is applicable to themethod being executed, the result of the system or method which isoutput can be at least one value of the one or more hardeningparameters, the value of the critical plastic strain, the value of atleast one parameter of the production function, or the generated sandproduction predictions.

Systems and methods also are provided for sand production from ageomechanical reservoir system, comprising operating the geomechanicalreservoir system in accordance with a result of any of the methods andsystems disclosed herein, including the foregoing systems and methods.

Systems and methods also are provided for operating a geomechanicalreservoir system to control sand production from the geomechanicalreservoir system. The systems and methods comprise computing a value ofat least one operating parameter based on a fit of a geomechanical modelcomprising a production function to data indicative of physical orchemical properties associated with materials within the reservoirsystem; wherein the production function predicts sand production fromthe geomechanical reservoir system; and wherein the at least one valueof the operating parameter indicates a condition for stabilized sandproduction from the geomechanical reservoir system; and operating thegeomechanical reservoir system in accordance with the value of the atleast one operating parameter.

4. BRIEF DESCRIPTION OF THE DRAWINGS

FIG. 1 is a flow chart of a method for use in predicting sand productionfrom a geomechanical reservoir system.

FIG. 2A depicts stress components (σ_(x), σ_(y), σ_(x)) for deviatoricloading of a sample 200 in a triaxial test (in this example,σ_(x)=σ_(z)).

FIG. 2B depicts the stress component (σ_(y)) for an oedometer test, inwhich a sample 202 is confined radially in a confining ring 204, and isloaded only in the y-direction (σ_(y)).

FIG. 3 depicts an example schematic of a hollow cylinder test.

FIG. 4 shows the results of a hollow cylinder test, plotted as sandproduction (in cubic centimeters (cc)) versus time (seconds) and a fitto the results using material constants of critical strain value=0.014,and the exponent (m)=1. The hollow cylinder sanding test data was fitfor both the massive sanding and the gradual sanding portions of thecurve.

FIG. 5A shows an example schematic of a perforation test.

FIG. 5B shows core samples from a perforation test.

FIG. 6 shows the results of a perforation test, plotted as sandproduction (in cubic centimeters (cc)) versus time (seconds) and fits tothe results using different values of the exponent (m). The cylindricalcore sample was perforated with a hemispherical end.

FIG. 7 is a flow chart of a method for use in operating a geomechanicalreservoir system based on a result of a sand production prediction.

FIG. 8 is a block diagram of an example approach for use in modelingsand production from a geomechanical reservoir system, including areservoir model and a geomechanical model.

FIG. 9 is a block diagram of an example approach for use in modelingsand production from a geomechanical reservoir system, including areservoir model, a geomechanical model, and a thermal model.

FIG. 10A illustrates an example of a three-dimensional grid for use incomputation of the models.

FIG. 10B illustrates an example of a two-dimensional grid for use incomputation of the models.

FIG. 11A shows a plot of the movement of the yield surface in aDrucker-Prager hardening model with shear hardening.

FIG. 11B shows a plot of the movement of the yield surface in aDrucker-Prager hardening model with cap hardening.

FIGS. 12A-D show plots of the yield surfaces in the octahedral(deviatoric) plane for the modified Drucker-Prager model for fourdifferent values of the deviation (K).

FIGS. 12E-H show plots of the yield surfaces in the octahedral(deviatoric) plane for the modified Bigoni-Piccolroaz model for fourdifferent combinations of values of beta (β) and the deviation (γ).

FIG. 13 illustrates an example computation of sand production using themodels.

FIG. 14 illustrates an example computer system for implementing the sandproduction prediction methods disclosed herein.

FIGS. 15A and 15B show core samples, including the dimensions of thesample, used for a brine test (using a hollow cylinder test setup suchas illustrated in FIG. 3).

FIG. 16A shows plots of the triaxial test data (black) and the numericalresults (data fit) from the simulation (grey) for the brine saturatedcore samples.

FIG. 16B shows a close-up of a portion of the plot shown in FIG. 16Awhere the core sample material becomes almost perfectly plastic.

FIG. 17 shows plots of the triaxial test results (black) and thenumerical results from the simulation (grey) for the kerosene saturatedsamples.

FIG. 18 illustrates the boundary conditions applied to the simulationsof the core samples to simulate the loading in a hollow cylinder test,where (A) shows the mechanical boundary conditions and (B) shows theflow boundary conditions, applied to the samples in the simulation.

FIG. 19 shows the measured and discretized confining stress, flow rateand flowing pressure vs. time for the brine saturated sample.

FIG. 20 shows a plot of the calculated permeability as a function offlow rate.

FIG. 21 shows a plot of the calculated permeability as a function ofconfining stress.

FIG. 22 shows a graph showing the match of the discretized and modeledflow rates as a function of time.

FIG. 23A shows the shape of the yield surface in the deviatoric planefor a Drucker-Prager-Type material model (beta=0 and deviation=0).

FIG. 23B shows the shape of the yield surface in the deviatoric planefor a Lade-Type material model (beta=0 and deviation=0.95).

FIG. 24 shows a plot of sand production volume (cc) vs. time (sec) froma measurement of a core sample (black line) and different fits of amodel to the results using grid sizes of 0.05 mm, 0.01 mm, and 0.005 mm.A grid size of 0.01 mm was small enough to capture the model details andeliminate grid dependencies. The material constants used were:deviation=0.5, beta=0, critical strain value=0.017, and the exponent=2.

FIG. 25 shows a plot of sand production volume (cc) vs. time (sec) froma measurement of a core sample (black line) and a fit to the resultsusing material constants of deviation=0.5, beta=0, critical strainvalue=0.017, and the exponent=2. The data was fit for both the massivesanding and the gradual sanding portions of the curve.

FIG. 26 shows a plot of the discretized measured water volume(liters(l)/min) vs. time (sec) from measurement of a core sample (blackline) and a fit to the results using material constants ofdeviation=0.5, beta=0, critical strain value=0.017, and the exponent=2.

FIG. 27 shows a plot of effective plastic strain vs. hardening (psi) forthe brine test.

FIG. 28 shows a plot of sand production volume (cc) vs. time (sec) frommeasurement of the core samples (black line) and different fits to theresults where the deviation and the critical strain limit are varied asindicated, and beta and the exponent are kept constant at 0 and 1,respectively.

FIG. 29 shows a plot of sand production volume (cc) vs. time (sec) frommeasurement of the core samples (black line) and different fits to theresults where the deviation=0.5, beta=0, and the critical strainvalue=0.017, and the exponent (m) value is varied as indicated.

FIG. 30 shows a plot of sand production volume (cc) vs. time (sec) frommeasurement of the core samples (black line) and different fits to theresults where the deviation=0.5, beta=0, and the exponent=1, and thecritical strain limit is varied.

FIG. 31 shows horizontal slices taken from the top to the bottom of acore sample, showing cavity progression. The slices were taken at theend of a sand production test. It was

FIG. 32 shows numerical simulations of the cavity progression over time.Each plot is a axially-symmetric vertical slice, where the grey areasshow the location and extent of the cavity and the black areas showintact rock (or spacer). Material constants used in the simulation were:deviation=0.5, beta=0, critical plastic strain value=0.017, and theexponent=2.

FIG. 33A shows a plot of sand production (cc) vs. time (sec), computedusing material constants of deviation=0.5, beta=0, critical strainvalue=0.017, and the exponent=2.

FIG. 33B shows a numerical simulation of cavity size, computed for aconstant confining stress of 46 MPa and a constant pressure differentialof 0.1643398 MPa from 116675 seconds onwards, and using materialconstants of deviation=0.5, beta=0, critical strain value=0.017, and theexponent=2. The grey area indicates the location of the cavity and theblack region indicates intact rock (or spacer).

FIG. 34 shows a numerical simulation of sand production (cc) vs. time(sec), computed for the a constant confining stress of 44 MPa and aconstant pressure differential of 1.255317 MPa from 113230 secondsonwards, and using material constants of deviation=0.5, beta=0, criticalstrain value=0.017, and the exponent=2.

5. DETAILED DESCRIPTION

Systems and methods are provided for use in predicting sand productionfrom a geomechanical reservoir system. The systems and methods usemeasured reservoir properties to generate such sand productionpredictions, which are useful in that they would allow determining,early in the reservoir development project, if and when sand productioncould occur and at what rate.

Such sand production predictions can be used to determine, early in areservoir development project, the type of completion techniques whichmay be implemented to drastically reduce sand production throughout thelife of a well, or to allow a certain amount of gradual sanding (agradual temporal evolution of the sand production), but which does notlead to massive sanding. For example, decisions can be made as to thedesigns of the barrier system to be used around the borehole of thereservoir, such as but not limited to sand screen technology or gravelpacks. Decisions can be made as to whether an unconsolidated formationwith multiple producing zones should be completed cased hole. Such sandproduction predictions could help to reduce the costs associated withwell completion, since the decision can be made, depending on theallowable degree of sand production, to install less sand productionmitigation equipment than would be used in the absence of the sandproduction predictions. The sand production predictions also can be usedto make decisions on how to operate the reservoir, such as but notlimited to the drawdown pressure, production rate, minimum bottomholepressure, temperature of the production zone, and fluid flow pressure inthe wellbore, to achieve the desired amount of sand production. Inaddition, the sand production predictions can be used to make a decisionas to the point in the lifetime of a well at which to use techniques andto install equipments to mitigate sand production.

The sand production prediction systems and methods disclosed herein arephysically based and can be used to predict the phenomena of onset ofsanding, gradual sanding, and massive sanding. The sand productionprediction systems and methods are deterministic, i.e., they are basedon physical principles and not merely correlations. For example,physically based sanding criteria are applied to the computations todetermine whether sanding occurs, in what amount, and at what rate.

Operating a reservoir according to the results of the sand productionpredictions can lead to a substantial savings by reducing the completioncosts. The results from the sand production predictions can be used forsanding management. Sanding management can be used, for example, to makedecisions concerning the type and cost of sand prevention completions toinstall, e.g., cased hole and perforation, or conventional proppedfracture completions, which produce a limited amount of sand.

5.1 Examples of Sand Production Prediction Systems and Methods

The flow chart of FIG. 1 shows steps in an example system and method foruse in predicting sand production from a geomechanical reservoir system.

Step 100. In step 100, data indicative of physical or chemicalproperties of material within the geomechanical reservoir system arereceived. For example, data indicative of physical or chemicalproperties of materials in the geomechanical reservoir system during oneor more stages of sanding of material may be received. In the example ofFIG. 1, data indicative of gradual sanding of material in thegeomechanical reservoir system are received in step 100A. Gradualsanding refers to a gradual temporal evolution of sanding of material inthe reservoir. In step 100B, data indicative of massive sanding ofmaterial in the geomechanical reservoir system are received. Massivesanding refers to an sharp increase in sand production. Such datainclude but are not limited to a rate of sanding (in units of volume ofsanding over time). Other examples of such data include, but are notlimited to, type of formation material, porosity of formation material,permeability of formation material, types of fluid in the reservoir,pore pressure, temperature, and viscosity of a fluid in the wellbore,temperature of the production zone, fluid flow pressure in the wellbore,drag force of fluid flow in the wellbore, and type of fluid flow in thewellbore. Other examples of such data include, but are not limited to,the depth of the well, the oil gravity, the oil viscosity, the netthickness of the rock type, the current reservoir pressure, the minimumoil content, the oil saturation, the permeability and porosity of therock, the temperature of the system, the transmissibility of the rockformation, the water salinity, existing fracture system, gas cap, dipangle, well spacing, receptivity, hydrocarbon (HC) composition, minimummiscibility pressure, pressure ratio, gas saturation, bubble pointpressure, critical gas saturation, gas ratio, and vertical sweep factor.

The received data may be indicative of plastic deformation of materialin the geomechanical reservoir system. That is, the data can be one ormore values of parameters which provide a measure of the plasticdeformation of formation materials. Data indicative of elastic behaviorof materials in the reservoir system also may be received. Examples ofsuch data are, but are not limited to, Young's modulus, yield strength,a stress-strain curve for the material, an ultimate strength, strainhardening behavior, necking behavior, point of fracture. In an example,data indicative of the onset of plastic deformation (that is, the pointof transition from elastic to plastic behavior) can be received. Forexample, yield strength can be used to pinpoint an elastic limit of amaterial, beyond which additional stress on the material can causepermanent (plastic) deformation to occur.

Determination of a yield criterion of a material also can be used todetermine an onset of plastic deformation of the material. A yieldcriterion (which can be displayed as a yield surface) can be used toindicate a limit of material elasticity (and an onset of plasticdeformation) under different combinations of stresses. Examples of yieldcriteria which can be as applied to isotropic materials, i.e., materialswhich have uniform properties in all directions, are criterion based ona maximum principal stress, a maximum principal strain, a maximum shearstress, a total strain energy, and a distortion energy. For a yieldcriterion based on a maximum principal stress, yield can be consideredto occur when the largest principal stress applied to the materialexceeds the uniaxial tensile yield strength. With a maximum principalstrain criterion, yield can be considered to occur when the maximumprincipal strain of the material reaches the strain corresponding to theyield point during a simple tensile test. For a maximum shear stressyield criterion (Tresca yield criterion), yield can be considered tooccur when the shear stress applied to the material exceeds the shearyield strength. In a total strain energy yield criterion, it can beassumed that the stored energy associated with elastic deformation atthe point of yield is independent of the specific stress tensor, so thatyield occurs when the strain energy per unit volume is greater than thestrain energy at the elastic limit in simple tension. With thedistortion energy yield criterion (Von Mises yield criterion), yield canbe considered to occur when the shape distortion of the material exceedsthe yield point for a tensile test. Other examples of yield criteriaapplied to isotropic materials are the Mohr-Coulomb yield criterion, theDrucker-Prager yield criterion, and the Bresler-Pister yield criterion.Examples of yield criteria which can be as applied to anisotropicmaterials, i.e., materials whose plastic yield behavior showsdirectional dependency, include but are not limited to, Hill's quadraticyield criterion, generalized Hill yield criterion, and the Hosford yieldcriterion.

Data indicative of elastic behavior, plastic deformation behavior,gradual sanding, or massive sanding, of a material in the reservoirsystem can be obtained from tests such as, but not limited to, atriaxial compression test, a triaxial extension test, a uniaxial straintest, an oedometer test, and a hydrostatic compression test. FIG. 2Ashows the direction and relationship of the stress components (σ_(x),σ_(y), σ_(z)=σ_(x)) which can be applied to a cylinder of material 200in a triaxial test. FIG. 2B shows the stress component (σ_(y)) which canbe applied in the y-direction to a cylinder of material 202 in anoedometer test; the material can be confined in a confining ring 204during the compression loading such that no stress is applied in thex,z-directions. Either or both of these tests can be used to providedata indicative of elastic behavior, plastic deformation behavior,gradual sanding, and/or massive sanding, of a material in the reservoirsystem. The data can be obtained from tests performed on the materialinvolving different loading paths for the material. Data also can beobtained from a hollow cylinder test, in which a flow test and differentcombinations of axial, spherical, or torsional stress are applied to ahollow cylindrical sample of material. FIG. 3 shows a schematic of ahollow cylinder test setup. FIG. 4 shows data obtained from a hollowcylinder test performed on a material, and also shows the point of onsetof plastic deformation, the region of gradual sanding, and the region ofmassive sanding in the data. Data also can be obtained from aperforation test, in which a flow test is performed on a cylindricalsample of a material after one end of it has been perforated (i.e.,subjected to a shaped explosive charge). Flow tests can be performedwith oil, gas, brine water, or any combination these fluids. FIG. 5Ashows a schematic of a perforation test setup, and FIG. 5B shows samplesof materials which have been subjected to a perforation test. FIG. 6shows data obtained from a hollow cylinder test performed on a material,and also shows the point of onset of plastic deformation, the region ofgradual sanding, and the region of massive sanding in the data.

Data can be obtained from tests performed on one or more core samplestaken from a site of an actual well, or the site of a proposed well. Forexample, data can be obtained from tests performed on a core sample ofmaterial can be taken from a wellbore. Such core samples can be takenfrom different parts of the reservoir representative of the formationwhere sanding may occur.

In another example, data can be obtained from tests performed on asynthetic sample which is created to have similar physical and chemicalproperties to the actual formation materials from a well, such as fromthe regions of the well formation where sanding may occur.

Step 102. Values of one or more hardening parameters are derived basedon a fit of at least one hardening model to the received data. Hardeningmodels can be used to model the plastic deformation behavior of materialwithin the geomechanical reservoir system. Hardening models arediscussed in Section 5.5.2.4 below. Examples of hardening models includethe Drucker-Prager model with shear hardening or cap hardening, themodified Drucker-Prager model with tabular hardening, the modifiedBigoni-Piccolroaz model, and the Matsuoka-Nakai model. Additionalexamples of hardening models include, but are not limited to, themodified Cam-Clay model, the DiMaggio and Sandler generalized cap model,the Lade model, the Iwan/Mroz multi-surface model, and the Fossum andFredrich continuous surface cap plasticity model.

Examples of hardening parameters which may be derived in connection withthe modified Drucker-Prager model are, but are not limited to, α (amultiplier of a first invariant of total stress), a Drucker-Pragerexponent, a yield constant (Γ), an effective plastic strain ε^(p), and adeviation (K). Examples of hardening parameters which may be derived inconnection with the modified Bigoni-Piccolroaz model are, but are notlimited to, deviation (γ), beta (β), α, and a yield constant (Γ).

Step 102 can involve several steps, as illustrated in FIG. 1. In step102A, a hardening model is selected. In a next step, the selectedhardening model is fit to the received data. For example, as illustratedin step 102B, the selected hardening model can be fit to received dataindicative of the plastic deformation behavior of the material in thereservoir (discussed above in connection with step 100). In step 102C,the values of one or more parameters of the hardening model fit to thedata are obtained.

The fit of the hardening model to the received plastic deformation datacan be performed using any applicable data fitting method. For example,the fit can be performed using a regression method, such as a linearregression and a nonlinear regression. Regression packages which performa regression fit to data are known in the art. The regression can beperformed with limited dependent variables, can be a Bayesian linearregression, a quantile regression, or a nonparametric regression. Thefit can be performed using a statistical method. Examples of packageswhich can be used for performing a fit of the hardening model to thereceived plastic deformation data include, but are not limited to, SAS®(SAS Institute Inc., Cary, N.C.), MATLAB® (The MathWorks, Inc., Natick,Mass.), R (accessible via the World Wide Web at the website of the RProject for Statistical Computing), and Dap (accessible via the WorldWide Web at the website of the GNU Operating System).

Step 104. In step 104, a value of a critical plastic strain is computedbased on a fit of a geomechanical model comprising a hardening model tothe received sanding data. In an example, the sanding data may be thereceived data indicative of gradual sanding data or the received dataindicative of massive sanding data, or a combination of the two. Thecritical plastic strain can be a critical value of the plastic strain ofthe material which indicates a point when the material fails, i.e.,rubblizes to form sand and generate a cavity. In an example, thecritical plastic strain is a critical value of an effective plasticstrain of the material.

An applicable geomechanical model can model the stresses, strains,and/or displacements of isotropic materials, transversely isotropicmaterials, linear elastic materials, porous materials, solid materials,or any combination thereof. In an example, the geomechanical model canmodel plastic deformation behavior of materials. Examples ofgeomechanical models are discussed in Section 5.5.2 below (whichincludes a discussion of hardening models).

The fit of the geomechanical model comprising the hardening model to thesanding data can be performed using any applicable data fitting method.For example, the fit can be performed using a regression method, such asa linear regression and a nonlinear regression. In another example, thefit can be performed by solving a system of partial differentialequations, where the system of partial differential equations comprisesthe geomechanical model comprising the hardening model. The system ofpartial differential equations also may comprise a reservoir flow model(discussed in Section 5.5.1 below). The system of partial differentialequations can be coupled through a fully-expanded Jacobian, as discussedin Section 5.7 below. The solving of the system of partial differentialequations can include solving simultaneously in a single time step thefully-expanded Jacobian based on the received data (see Section 5.7below), such as the sanding data. The system of partial differentialequations also may comprise a thermal model (discussed in Section 5.5.3below). A procedure which can be used for solving the system of partialdifferential equations is discussed below in Section 5.7.

The fit of the geomechanical model comprising the hardening model to thesanding data can be performed as part of a broader computation of asystem of equations which model the geomechanical reservoir system, andwhich can be used to compute the stresses, strains, and/or displacementsthat arise when fluids are injected into or produced from a reservoir aswell as when stresses are applied to the boundaries of a reservoir. Insome computations, a reservoir system model, comprising a geomechanicalmodel (comprising the hardening model), a thermal model, and a reservoirflow model, is capable of solving systems that include porous flow, heatflow, and geomechanics.

In an example, the fit of the geomechanical model comprising a hardeningmodel to the received massive sanding data can be obtained by solving asystem of partial differential equations that model the geomechanicalreservoir system, where the system of partial differential equationscomprises a reservoir flow model and the geomechanical model comprisingthe hardening model, where the system of partial differential equationsis coupled through a fully-expanded Jacobian, and where solving thesystem of partial differential equations includes solving simultaneouslyin a single time step the fully-expanded Jacobian based on the receivedmassive sanding data. The system of partial differential equations alsomay comprise a thermal model.

Although step 102 was discussed prior to step 104, it should be notedthat either step could have been performed first. As illustrated in FIG.1, either of steps 102 and 104 may be performed at that stage of theflow chart. Steps 102 and 104 also may be performed iteratively, where aresult from performance of one step can be used in the other step. Forexample, if step 102 is performed first, then the value of one or moreparameters from the fit in step 102 is used in step 104 to compute thecritical plastic strain. If step 104 is performed first, then thecritical plastic strain computed in step 104 is used in the fit in step102 to derive values of one or more hardening parameters. Furthermore,steps 102 and 104 may be performed repeatedly and iteratively to arriveat better fits to the data. As a non-limiting example, if the hardeningmodel is a modified Bigoni-Piccolroaz model, then an initial value ofone or more of the deviation (γ), beta (β), α, or yield constant (Γ) canbe derived in step 100 and used in step 104 to compute, for example, avalue of a critical plastic strain which fits to data indicative ofmassive sanding. The computed value of the critical plastic strain maybe used in a second iteration of computation in step 100 of a parameterof the hardening model. The second iteration of Step 100 could provide abetter fit to the received data. The result from the second iteration ofstep 100 may be provided to step 104 for a second iteration ofcomputation of the critical plastic strain to improve the fit to themassive sanding data. The iterations may cease once the fits in steps100 and 104 have converged to a best match to the data.

Step 106. A value of at least one parameter of a production function iscomputed based on a fit of a geomechanical model comprising theproduction function to the received gradual sanding data, and using thevalue of the critical plastic strain (discussed in step 104) and a valueof at least one of the hardening parameters (discussed in step 102).

The production function models the amount of sand production from amaterial prior to reaching the critical plastic strain. Productionfunctions are discussed in Section 5.5.2.5 below. The productionfunction can be any function ƒ(x) whose values vary from ƒ(x)=0 whenx=0, to ƒ(x)=1 when x=1. The sand production function also can be anyfunction ƒ(x) whose values can be scaled to fall within a range fromƒ(x)=0 when x=0 and ƒ(x)=1 when x=1. In another example, the sandproduction function also can be any function ƒ(x) whose values can betransformed to fall within the range from ƒ(x)=0 when x=0 and ƒ(x)=1when x=1, by application of a suitable transform, such as but notlimited to, a wavelet transform and a Lagrange transform. The term x canbe any function g( ) of the critical plastic strain (i.e., x=g(ε_(lim)^(p))). The function ƒ(x) may be any monotonic function of x, includingbut not limited to, a fraction function, a power function, a sinefunction, a cosine function, a logarithmic function, an exponentialfunction, a sigmoid function, or any combination thereof. In someexamples, x can be a function of both the critical plastic strain(ε_(lim) ^(p)) and the plastic strain invariant (ε^(p)) of a material,such as but not limited to a ratio x=ε^(p)/ε_(lim) ^(p). In an example,the production function can be given by ƒ(x)=(ε^(p)/ε_(lim) ^(p))^(m),where m is an exponent.

The parameter for which a value is computed in step 106 may be anyparameter which can be used to characterize the production function. Forexample, if the production function is a power function, the parametercan be at least one exponent of the power function. In a certain examplewhere the production function is given by ƒ(x)=(ε^(p)/ε_(lim) ^(p))^(m),the parameter can be the exponent m. In the foregoing example, the atleast one parameter of the production function for which a value is acomputed in step 106 is an exponent. In other examples, the parametercan be a multiplier of a term of the production function.

The production function can be used to predict the amount of sandproduction from a material prior to failure, i.e., during gradualsanding prior to reaching the critical plastic strain.

For example, FIG. 4 shows the results of a fit of the productionfunction (ε^(p)/ε_(lim) ^(p))^(m) to the data received from a hollowcylinder test. The data was fit for both the massive sanding and thegradual sanding portions of the curve using material constants ofcritical strain value (ε_(lim) ^(p))=0.014, and the exponent (m)=1. Inanother example, FIG. 6 shows the results of different fits of theproduction function (ε^(p)/ε_(lim) ^(p))^(m) to the data received from aperforation test using different values of the exponent (m).

The geomechanical model comprising the production function can be fit tothe sanding data using any applicable data fitting method. For example,the fit can be performed using a regression method, such as a linearregression and a nonlinear regression. In another example, the fit canbe performed by solving a system of partial differential equations,where the system of partial differential equations comprises thegeomechanical model comprising the production function. The system ofpartial differential equations also may comprise a reservoir flow model.The system of partial differential equations can be coupled through afully-expanded Jacobian and the solving of the system of partialdifferential equations can include solving simultaneously in a singletime step the fully-expanded Jacobian based on the received data, suchas the sanding data. The system of partial differential equations alsomay comprise a thermal model.

In an example, the fit of the geomechanical model comprising theproduction function to the received gradual sanding data can be obtainedby solving a system of partial differential equations that model thegeomechanical reservoir system, where the system of partial differentialequations comprises a reservoir flow model and the geomechanical modelcomprising the production function, where the system of partialdifferential equations is coupled through a fully-expanded Jacobian, andwhere solving the system of partial differential equations includessolving simultaneously in a single time step the fully-expanded Jacobianbased on the received gradual sanding data.

The fit of the geomechanical model comprising the production function tothe sanding data can be performed as part of a broader computation of asystem of equations which model the geomechanical reservoir system, andwhich can be used to compute the stresses, strains, and/or displacementsthat arise when fluids are injected into or produced from a reservoir aswell as when stresses are applied to the boundaries of a reservoir. Insome computations, a reservoir system model, comprising a geomechanicalmodel (comprising the production function), a thermal model, and areservoir flow model, is capable of solving systems that include porousflow, heat flow, and geomechanics.

In step 108, information which can be used to predict sand productionfrom the geomechanical reservoir system may be output. Such informationcan be, but is not limited to, at least one value of the one or morehardening parameters, the value of the critical plastic strain, and/orthe value of at least one parameter of the production function. Theinformation can be output to a user, or to various components, such asto a user interface device, a computer readable storage medium, amonitor, a user-accessible local computer, or a user-accessible computerthat is part of a network. For example, the output can be visuallydisplayed to a user using a monitor or user interface device such as ahandheld graphic user interface (GUI) including a personal digitalassistant (PDA).

5.1.1 Other Sand Production Prediction Systems and Methods

Other examples systems and methods for use in predicting sand productionfrom a geomechanical reservoir system includes a step of applying a sandproduction criterion to a geomechanical model which comprises one ormore hardening models. The system and method comprises the steps ofreceiving data indicative of physical or chemical properties associatedwith materials within the geomechanical reservoir system, and generatingsand production predictions by solving a system of partial differentialequations that model the geomechanical reservoir system. In addition tothe geomechanical model, the system of partial differential equationsmay comprise a reservoir flow model and/or a thermal model of thegeomechanical reservoir system.

One or more sand production criteria may be applied to the geomechanicalmodel. The sand production criteria can be determined as when a criticalvalue is reached for: (1) the total strain invariant, or (2) the plasticstrain invariant, or (3) the maximum effective stress. When the criticalvalue of the sand production criterion is reached, the material fails,i.e., rubblizes to form sand and generate a cavity. The sand productioncriteria are discussed in Section 5.5.2.5 below.

In an example, the system of partial differential equations can becoupled through a fully-expanded Jacobian, where the solving of thesystem of partial differential equations, such as by using a computersystem, includes solving simultaneously in a single time step thefully-expanded Jacobian based on the received data. The generated sandproduction predictions may be output to a user, a user interface device,a monitor, a computer readable storage medium, a local computer, or acomputer that is part of a network. For example, the generated sandproduction predictions can be visually displayed to a user using amonitor or user interface device such as a handheld graphic userinterface (GUI) including a personal digital assistant (PDA).

5.2 Systems and Methods For Operating a Geomechanical Reservoir System

Systems and methods for use in controlling sand production from ageomechanical reservoir system during operation also are disclosed. Themethod operating the geomechanical reservoir system in accordance with aresult of the implementation of any of the sand production predictionsystems and methods disclosed herein. The flow chart of FIG. 7 showssteps in an example system and method for use in operating ageomechanical reservoir system based on a result of a sand productionprediction.

In step 700, data indicative of physical or chemical propertiesassociated with materials within the reservoir system are received. Thedata received in step 700 can include any of the data described in step100 above.

In step 702, a value of at least one operating parameter is computedbased on a fit of a geomechanical model comprising a production functionto the received data, where the production function predicts sandproduction from the geomechanical reservoir system, and where the atleast one value of the operating parameter indicates a condition forsand production from the geomechanical reservoir system. In an example,the computed value of the operating parameter indicates a condition forstabilized sand production from the geomechanical reservoir system.

The prediction of sand production from a reservoir from any of thesystems or methods disclosed herein may be used to compute or derive avalue of at least one operating parameter for operating a reservoir toachieve the desired amount of sand production or the desired sandproduction behavior. Examples of operating parameters include but arenot limited to, the drawdown pressure, production rate, minimumbottomhole pressure, temperature of the production zone, fluid flowpressure in the wellbore, confining stress, and pressure differential.The results of implementation of the sand production predictions may beused to compute the value of one or more operating parameters whichcause a substantially stabilized cavity to develop and maintain in thereservoir. A cavity in a reservoir can be substantially stabilized if,over time, the cavity does not grow substantially larger, or growslarger to a negligible extent. For example, values of one or moreoperating parameters may be derived based on results of the sandproduction predictions which indicate an operating condition for thereservoir whereby a limited amount of sand is produced from thereservoir and a cavity in the reservoir created by that limited amountof sanding is substantially stabilized. The sand production from areservoir may be kept to a minimum by controlling production variablessuch as, but not limited to, drawdown, minimum bottomhole pressure, andproduction rate.

In step 704, the geomechanical reservoir system is operated inaccordance with the value of the at least one operating parameter. Thatis, the geomechanical reservoir system can be operated at the computedvalue for the at least one operating parameter, such as a value ofdrawdown pressure, production rate, minimum bottomhole pressure,temperature of the production zone, fluid flow pressure in the wellbore,confining stress, pressure differential, or any combination of theseparameters.

The results from the implementation of any of the sand productionprediction systems and methods can be used to determine the type ofcompletion techniques which may be installed in a reservoir to achievethe desired amount of sand production throughout the life of a well. Forexample, the designs of the barrier system to be used around theborehole of the reservoir (such as but not limited to sand screentechnology or gravel packs) may be selected based on these results. Thepredictions of sand production also could influence the well completionstrategy, such as open-hole, cased hole, and perforated cased hole, useof screen liners, or frac-and-pack (hydraulic fracturing followed byinjection of a proppant into the fracture). In addition, the sandproduction predictions can be used to make a decision as to the point inthe lifetime of a well at which to use techniques and to installequipments to mitigate sand production.

5.3 Examples of Modeling Methods

FIGS. 8 and 9 depict examples of systems for use in predicting sandproduction from a geomechanical reservoir system during oil production.Applicable modeling system includes one or more models which describevarious physical aspects of a geomechanical reservoir system. FIG. 8 andFIG. 9 depict modeling systems which include a reservoir fluid flowmodel and a geomechanical model. The reservoir fluid flow modeldescribes, e.g., porous flow, production and injection. Thegeomechanical reservoir model describes, e.g., stresses, strains, anddisplacement that arise when fluids are injected into or produced from areservoir and when stresses are applied to the boundaries of areservoir. The system of FIG. 9 also includes a thermal model. Thethermal model described heat flow. A system of non-linear partialdifferential equations can be used to interrelate the various aspects ofthese models.

After receiving data indicative of physical or chemical properties ofmaterial within the geomechanical reservoir system (such as but notlimited to plastic deformation, gradual sanding or massive sanding), asolver generates predictions (e.g., sand production predictions) byapplying the steps of the method in Section 5.1 above, and at pertinentpoints, by solving the system of partial differential equations. In thesolver of FIGS. 8 and 9, the system of partial differential equationscan be coupled through a fully-expanded Jacobian. The solving of thesystem of partial differential equations includes solving simultaneouslyin a single time step the fully-expanded Jacobian based upon thereceived data.

As depicted in FIGS. 8 and 9, the sand production prediction steps(steps 100 to 106 discussed above) can be performed iteratively with thesolution of the non-linear system of partial differential equations.That is, at the steps of the sand production prediction which includingperforming a fit, for example, each of steps 102, 104 and 106, thecomputation can be performed by solving simultaneously in a single timestep a fully-expanded Jacobian of the equations representing thegeomechanical reservoir system which are pertinent to the particularstep. The generated sand production predictions can be output to variouscomponents, such as output to a user interface device, a computerreadable storage medium, a monitor, a user-accessible local computer, ora user-accessible computer that is part of a network.

The non-linear system of partial differential equations compriseequations that correspond to what models are to be used in analyzing thegeomechanical reservoir system in the given step of the sand productionprediction. For example, FIG. 8 provides an example where the non-linearsystem of partial differential equations includes equationscorresponding to a reservoir flow model and a geomechanical model of thegeomechanical reservoir system. Depending on the step of the sandproduction prediction being performed, the geomechanical model mayinclude one or more hardening models, a sand production model(comprising the production function), or both. As another illustration,FIG. 9 provides an example where the non-linear system of partialdifferential equations includes equations corresponding to a reservoirflow model, a geomechanical model, and a thermal model of thegeomechanical reservoir system. Examples of equations that correspond toeach of the different models of the geomechanical reservoir system areprovided in Section 5.5.

Coupling of the various aspects of the models can be implemented such asthrough variables in a fully-expanded Jacobian. For example, afully-expanded Jacobian can act to couple fluid flow in the reservoir tothe geomechanical model by one or more of the following variables:effective stress, a porosity and one or more displacements associatedwith the geomechanical model. A variable in the fully-expanded Jacobianthat couples the geomechanical model to the fluid flow can be porosityand permeability that is associated with the reservoir flow model. Avariable in the fully-expanded Jacobian that couples the thermal modelto the geomechanical model can be a thermal stress associated with thethermal model. A variable in the fully-expanded Jacobian that couplesthe thermal model to the reservoir flow model can be fluid viscosity,conduction and convection in the reservoir associated with the thermalmodel. The fully-expanded Jacobian may include terms related to a rateof change (i.e., a time derivative), a spatial derivative, or partialspatial derivative, of a coupling variable, where the derivatives can beof any order, for example, a first-order derivative, a second-orderderivative, a third-order derivative, etc. First-, second-, third-and/or higher-order derivatives (whether time derivatives or spatialderivatives) of the coupling variables can be included in thefully-expanded systems of equations. Examples of variables that maycouple the different models are provided in Section 5.6.

The nonlinear system of equations of the fully-expanded Jacobian can besolved through numerical approaches, such as the approach discussed ingreater detail in Section 5.7, wherein the nonlinear system of equationsis solved, e.g., using a full Newton-Raphson expansion of all solutionvariables, which enhances solution stability and allows second orderconvergence rates for the nonlinear iterations. Examples of apparatusand computer-program implementations of the different methods disclosedherein are discussed in Section 5.8.

In another aspect, a system and method can include the steps ofreceiving data indicative of physical or chemical properties associatedwith material within the geomechanical reservoir system, defining a gridcomprising a plurality of grid cells, and generating sand productionpredictions by solving a system of partial differential equations thatmodel the geomechanical reservoir system. In addition to thegeomechanical model, the system of partial differential equations maycomprise a reservoir flow model and/or a thermal model of thegeomechanical reservoir system.

One or more sand production criteria may be applied to the geomechanicalmodel. The sand production criteria can be determined as when a criticalvalue is reached for: (1) the total strain invariant, or (2) the plasticstrain invariant, or (3) the maximum effective stress. The sandproduction criteria are discussed in Section 5.5.2.5.

In computations which use criterion (2), i.e., computations involvingcomputation of the critical value of the plastic strain invariant (thecritical plastic strain), steps 100-106 may be implemented, and at leastone parameter of a production function may be computed. The productionfunction can be used to predict the amount of sand production from agrid cell prior to reaching the critical plastic strain value.

In an example, the system of partial differential equations can becoupled through a fully-expanded Jacobian, where the solving of thesystem of partial differential equations, such as by using a computersystem, includes solving simultaneously in a single time step thefully-expanded Jacobian based on the received data. The reservoir model,thermal model and the geomechanical model may be computed onthree-dimensional grid cells or two-dimensional grid cells. Three andtwo dimensional grid cells which may be used in the simulation methodsherein are described in Section 5.4.

The generated sand production predictions may be output to a user, auser interface device, a monitor, a computer readable storage medium, alocal computer, or a computer that is part of a network.

Solving the nonlinear system of equations implicitly, e.g., using a fullNewton-Raphson expansion of solution variables, can enhance numericalstability (e.g., when dealing with cavity generation, or with anysimulation that involves very small grid blocks). Using a fully expandedJacobian from an implicitly coupled system of equations provides morestability for the solution process.

5.4 Simulation Method

FIG. 10A illustrates an example of a three-dimensional (3D) grid whichcan be used for the computations of the geomechanical model, and thereservoir model and/or the thermal model. For example, one or more ofthe multi-point flux (such as for reservoir and porous flow) model, thegeomechanical model (such as for computing stress, displacement, and/orcavity generation (rubblizing)), and the thermal model, may be computedon the 3D grid. The computation of the geomechanical model, and thereservoir model and/or the thermal model can use a parallel processingapproach to couple the grids. The dimensions of the grid cells can be onthe order of feet, inches, or fractions of an inch. In another example,the dimensions of the grid cells can be on the order of meters,centimeters, millimeters, microns, or fractions of a micron.

The 3D grid may be structured or unstructured hexahedral gridscomprising hexahedral elements. A hexahedral grid cell has eightcorners, twelve edges (or sides), and six faces. The hexahedral gridcells may each include at least eight nodes (one at each corner), ormore and up to twenty-seven (27) nodes (i.e., a node at the center ofeach face, the center of each side, the center of each edge, and in thecenter of the cell). Different hexahedral cells may include differentnumbers of nodes. In an example, the 3D grid may include structured orunstructured tetrahedral grid elements. In another example, the 3D gridmay include other morphologies of elements which span the range betweentetrahedral and hexahedral grid elements, either structured orunstructured. The 3D grid may include any combination of theaforementioned grid elements.

A two-dimensional (2D) grid also can be used for the computations of thegeomechanical model, and the reservoir model and/or the thermal model.For example, a 2D grid can be used for axisymmetric computations. FIG.10B illustrates an example of a two-dimensional (2D) grid which can beused for the computations of the reservoir model, thermal model and thegeomechanical model.

The 2D grid may be structured or unstructured quadrilaterals gridscomprising quadrilateral elements. Each quadrilateral grid cell has fourcorners and four edges. The quadrilateral grid cells each include atleast four nodes (one at each corner) and up to five nodes (i.e., a nodeat the center).

In certain examples, computations can be performed on both a 2D grid anda 3D grid. For example, depending on the topography of the reservoir,certain of the computations can be performed on a 2D grid, while othercomputations can be performed on a 3D grid. In these computations, thenodes of the 2-D grid can be configured to coincide with the nodes onone of the outer boundaries of the 3D grid and certain computations,such as fracture widths and 3D displacements, can be coupled at commonnode points. The input data format for the 2D grid can be similar tothat for the 3D grid.

For certain computations, parameters such as fluid flow, displacements,cavity generation, and tractions, can be monitored across the elementsof the grid cells.

5.5 Models

Examples of the differential equations that correspond to each of thedifferent models of the geomechanical reservoir system are providedbelow. The differential equations for the models included in acomputation can be combined to produce an implicit, fully-coupledformulation. A consistent set of units can be used for all variables inequations included in a computation.

5.5.1 Reservoir Model

The system of equations for porous flow include conservation of mass

$\begin{matrix}{{{\frac{\partial}{\partial t}({\rho\varphi})} = {{- {\nabla\left( {\rho \; \overset{->}{U}} \right)}} + q_{w}}},} & (1)\end{matrix}$

where φ is porosity and ρ is fluid density which may be a function ofpressure. The model allows wells to be completed in the reservoirelements and the q_(w) in the equation above accounts for injection intothe reservoir elements.

The velocity {right arrow over (U)} is the Darcy velocity relative tothe porous material and can be defined by

$\begin{matrix}{{\overset{\rightharpoonup}{U} = {{- \frac{K}{\mu}}\left( {{\nabla p} - {\rho \; g{\nabla h}}} \right)}},} & (2)\end{matrix}$

where K is a tensor permeability, μ is the viscosity which may befunction of pressure, p is fluid pressure, and ρg∇h is a gravitationalterm.

The geomechanical variables included in the fluid flow equationshighlight the coupling between the flow and deformation models (thedefinitions of some geomechanical terms are described in Section 5.5.2below).

Temperature-dependent water properties may be entered in different waysfor computations involving temperature changes. The water properties maybe entered as functions of both pressure (P) and temperature (T) for thefully coupled computations. For iteratively coupled computations, thewater properties can be entered as functions of pressure and thenmodification factors can be used for the temperature effects. Thetreatment for temperature dependent properties is explained in moredetail in Section 5.5.3.

The thermal behavior of the fluids also can be modeled by modifying thefluid properties using modification factors (described in Section5.5.3.3 below).

5.5.1.1 Multiphase Porous Flow

The reservoir model allows for several phase behavior models rangingfrom single phase to black-oil to fugacity-based compositions. Darcyflow may be modeled for aqueous phases, nonaqueous liquid phases, andnonaqueous vapor phases and nc components. Any phase behavior models maybe used with the porous flow models disclosed herein. The fluid flowequations may be presented in terms of a general compositionalformulation. Partial differential equations representing component massbalances for multiphase flow are:

$\begin{matrix}{{{\frac{\partial}{\partial t}\left( {\varphi N}_{ic} \right)} = {{- {\sum\limits_{\alpha}{{\nabla x_{ic}^{\alpha}}\rho_{\alpha}{\overset{\rightarrow}{U}}_{\alpha}}}} + q_{ic}}},} & (3)\end{matrix}$

where N_(ic) is the concentration of component ic per unit pore volume,given by

${N_{ic} = {\sum\limits_{\alpha}{x_{ic}^{\alpha}\rho_{\alpha}S_{\alpha}}}},$

x_(ic) ^(α) is the mole fraction of component ic in phase α, ρ_(α) isthe molar density of phase, and q_(ic) is the molar flow rate ofcomponent ic per unit reservoir volume. The velocity of phase α is givenby

$\begin{matrix}{{U_{\alpha} = {{- \frac{K\; k_{r\; \alpha}}{\mu_{\alpha}}} \cdot \left( {{\nabla P_{\alpha}} - {\gamma_{\alpha}{\nabla\; D}}} \right)}},} & (4)\end{matrix}$

Phase pressures can be defined by

P _(α) =P+P _(cα),  (5)

where P_(cα) is the capillary pressure and P is the reference pressure.The reference pressure is used for PVT calculations, well calculations,and geomechanical calculations. The reference pressure can be thenonaqueous phase pressure for two-phase models, and the nonaqueousliquid phase pressure for three-phase models.

The porosity can be defined as

φ=φ_(o)[1+c _(r)(P−P _(o))]  (6)

for porous flow models where φ_(o), c_(r), and the initial pressureP_(o) are functions of location.

For computations that include the geomechanical model, the porosityrelative to the initial undeformed bulk volume is given by Equation 15,where it is seen that c_(r) can be related to the Biot constant forporosity.

5.5.1.2 Non-Darcy Flow

In certain examples, non-Darcy flow may be modeled using theForschheimer equation to modify the relationship between the pressuregradient and the fluid velocity. In other example, non-Darcy flow can bemodeled by specifying a general relationship between fluid velocity andpressure gradient.

The Forschheimer equation for the non-Darcy velocity (which wouldreplace Equation 2 above), for systems that involve 3-D, multiphase flowin anisotropic media, can be given by:

$\begin{matrix}{{{{- \frac{k_{r\; \alpha}}{\mu_{\alpha}}}\underset{\_}{K}{\overset{\rightarrow}{\nabla}p_{\alpha}}} = {\left( {1 + {\beta_{K}\rho_{\alpha}\frac{k_{r\; \alpha}}{\mu_{\alpha}}{{\overset{\rightarrow}{v}}_{\alpha}}_{2}}} \right){\overset{\rightarrow}{v}}_{\alpha}}},} & (7)\end{matrix}$

where μ_(α) is the viscosity of phase α, K is the permeability tensork_(rα) is the relative permeability of phase α, v_(α) is the Darcyvelocity of phase α, the expression ∥ . . . ∥₂ is the L₂-norm defined by∥{right arrow over (u)}∥₂=√{square root over (u_(x) ²+u_(y) ²+u_(z) ²)},and the parameter β_(k) is a non-Darcy coefficient. Parameter β_(k) mayvary throughout the reservoir, therefore, a value for β_(k) may beentered for each grid block. The non-Darcy coefficient β_(k) is relatedto the inverse of the transition constant, i.e. β_(k)=1/Γ. See Barree etal., “Beyond Beta Factors, A Complete Model for Darcy, Forchheimer, andTrans-Forschheimer Flow in Porous Media,” SPE 89325, SPE annualTechnical Conference and Exhibition, Houston, Tex., Sep. 26-29, 2004.The Reynolds number for phase is given by the equation

$\begin{matrix}{N_{re}^{\alpha} = {\beta_{k}\rho_{\alpha}\frac{k_{r\; \alpha}}{\mu_{\alpha}}{{{\overset{\rightarrow}{v}}_{\alpha}}_{2}.}}} & (8)\end{matrix}$

The units of the terms in Equation 8 should be chosen such that theresult is dimensionless. After combining Equations 7 and 8, theForschheimer equation becomes

$\begin{matrix}{{{\overset{\rightarrow}{v}}_{\alpha} = {{- \frac{k_{r\; \alpha}}{\mu_{\alpha}}}\left( {1 + N_{re}^{\alpha}} \right)^{- 1}\underset{\_}{K}\; {\overset{\rightarrow}{\nabla}p_{\alpha}}}},} & (9)\end{matrix}$

where non-Darcy flow is expressed as a phase-dependent, permeabilitymodification function that varies with the Reynolds number for thatphase. This is different from the standard permeability modificationfunction, because this non-Darcy formulation has separate values foreach phase.

In some cases, the Forschheimer equation does not provide an adequateapproximation for non-Darcy flow. Modification functions of thefollowing form can be used to approximate non-Darcy flow:

$\begin{matrix}{{{\overset{\rightarrow}{v}}_{\alpha} = {{- \frac{k_{r\; \alpha}}{\mu_{\alpha}}}{f_{\alpha}\left( N_{re}^{\alpha} \right)}\underset{\_}{K}\; {\overset{\rightarrow}{\nabla}p_{\alpha}}}},} & \left( {9A} \right)\end{matrix}$

Modification functions can be constructed which satisfy the constraintsƒ_(α)(N_(re) ^(α))≦1.0 and ƒ_(α)(0)=1.0. For the standard Forschheimerequation, the following function can be specified:

ƒ_(α)(N _(re) ^(α))=(1+N _(re) ^(α))⁻¹,  (9B)

For an extended Forschheimer equation, the following functions can bespecified:

$\begin{matrix}{{{f_{\alpha}\left( N_{re}^{\alpha} \right)} = {K_{ratio} + \frac{1 - K_{ratio}}{1 + N_{re}^{\alpha}}}},} & \left( {9C} \right)\end{matrix}$

Equation 3 in J. L. Miskimins, et al., “Non-Darcy Flow in HydraulicFractures: Does It Really Matter?” SPE paper 96389, SPE annual TechnicalConference and Exhibition, Dallas, Tex., Oct. 9-12, 2005, is anotherform of a non-Darcy formulation which is applicable to the methodsdisclosed in this application.

5.5.1.3 Computation in the Reservoir Model

Computation of the reservoir model can be on the 3D grid (which can bethe grid used for the geomechanical model). The 3D reservoir grid caninclude porous flow calculations. The fluid velocity terms can becomputed for flow between reservoir cells as well as for flow betweenreservoir cells and cavity cells (cells where sanding has occurred). Aprimary variable for the porous flow equations can be the fluid pressureor the fluid composition, which can be evaluated at the center of eachhexahedral element (cell-based). In certain computations, a multi-pointflux algorithm can be used for the unstructured reservoir flowcomputations so the resulting computational stencil can be 27-point forgeneral hexahedral elements of a 3D grid when eight elements share acommon corner.

5.5.2 Geomechanical Model

5.5.2.1 Poroelastic Materials

A linear relationship (small strain) can be used for thestrain-displacement relations. The coupled flow/displacement modelrelating stress, strain, temperature, and pore pressures can be basedupon a Biot poroelastic theory. The equilibrium equation can be basedupon total stresses and assumes quasi-static equilibrium.

The poroelastic equations can be formulated in terms of total stresses,bulk strains, temperatures, and pore pressures. The total stress can bedefined by the average tractions that one would observe on a planarsection of the reservoir where the planar section includes loads carriedby the solid and pore pressures from the fluid. The bulk strains can bethe strains that one would observe from a strain gauge if it wereattached to the deforming porous material.

The system of equations for the linear poroelastic displacements includethe strain-displacement equation

$\begin{matrix}{{ɛ_{ij} = {\frac{1}{2}\left( {u_{i,j} + u_{j,i}} \right)}},} & (10)\end{matrix}$

where the commas imply differentiation, u_(i) is the displacement in thei-direction ε_(ij) is the bulk strain of the porous material, andexpansion corresponds to positive strains. The total stresses satisfythe equilibrium equations

σ_(ij,j)+ƒ_(i)=0  (11)

where the stress tensor is symmetric, and the gravity term f_(i) is afunction of the solid density, fluid densities, and porosity. Tractionor displacement boundary conditions can be specified in all threedirections at all six boundaries of the three-dimensional grid on whichthe model is computed.

When differences in temperature are not taken into account, theconstitutive equations relating total stresses, bulk strains, and porepressure are

$\begin{matrix}{{\sigma_{ij} = {\sigma_{ij}^{o} + {\frac{E}{1 + v}\left( {ɛ_{ij} + {\frac{v}{1 - {2v}}ɛ_{kk}\delta_{ij}}} \right)} - {{\alpha \left( {p - p_{o}} \right)}\delta_{ij}}}},} & (12)\end{matrix}$

where tension is positive, the repeated index kk implies summation,σ_(ij) ^(o) is the initial in-situ stress, p_(o) is the initialpressure, E is the elastic modulus, v is Poisson's ratio, α is Biot'sconstant in stress/strain equations, δ_(ij) is 1 when i=j and 0 wheni≠j. It can be assumed that the strains are zero when σ_(ij)=σ_(ij)^(o).

In the examples where differences in temperature are taken into account,the constitutive equations are:

$\begin{matrix}{{\sigma_{ij} = {\sigma_{ij}^{o} + {\frac{E}{1 + v}\left( {ɛ_{ij} + {\frac{v}{1 - {2v}}ɛ_{kk}\delta_{ij}}} \right)} - {{\alpha \left( {p - p_{o}} \right)}\delta_{ij}} - {\alpha_{T}{K\left( {T - T_{o}} \right)}\delta_{ij}}}},} & (13)\end{matrix}$

where α_(T) is the thermal volumetric expansion coefficient forstress/strain equations, and K is the elastic bulk modulus. The pressurep_(o) is the initial pore pressure and T_(o) is the initial temperature.

If the stress and pressure terms are combined to form σ_(ij)^(e)=σ_(ij)+αpδ_(ij), then the equation becomes a standard thermallinear elastic constitutive equation where the stresses have beenreplaced by the effective stresses σ_(ij) ^(e). If the initial in-situstresses and initial pore pressure are zero, the equation then takes thestandard form

$\begin{matrix}{{\sigma_{ij}^{e} = {{\frac{E}{1 + v}\left( {ɛ_{ij} + {\frac{v}{1 - {2v}}ɛ_{kk}\delta_{ij}}} \right)} - {\alpha_{T}{K\left( {T - T_{o}} \right)}\delta_{ij}}}},} & (14)\end{matrix}$

The relationship between the porosity (relative to the undeformed bulkvolume) and the strains and fluid pressure (when differences intemperature are not taken into account) is given by:

$\begin{matrix}{{\varphi = {\varphi_{o} + {\alpha \; ɛ_{ii}} + {\frac{1}{M}\left( {p - p_{o}} \right)}}},} & (15)\end{matrix}$

where Equation 15 assumes that the initial strains are zero, φ_(o) isthe initial porosity, and M⁻¹ is Biot's constant for pore pressure inporosity equations.

When differences in temperature and deposition fraction are taken intoaccount, the porosity relative to the undeformed bulk volume is definedas:

$\begin{matrix}{{\varphi = {\varphi_{o} + {\alpha \left( {ɛ_{v} - \; ɛ_{v}^{o}} \right)} + {\frac{1}{M}\left( {P - P_{o}} \right)} - {\alpha_{V}\left( {T - T_{o}} \right)} - \sigma}},} & (16)\end{matrix}$

where, α and M⁻¹ are Biot constants, P is the phase pressure (formultiphase flow), α_(V) is the thermal volumetric expansion coefficientfor porosity, and σ is the deposition fraction (the volume fraction ofsolid waste deposited per bulk volume, for example, of a grid element ina computation). Solid waste can be deposited within pores as waste movesthrough the reservoir, and there can be reductions in porosity andpermeability as the waste deposition builds up in pore spaces. Theporosity in a computation grid element can be a function of fluidpressure, temperature, and deformations, while the amount of porosityreduction due to deposition of solids can be set equal to σ.

For an isotropic material, the six poroelastic material parameters: E,v, α, and M⁻¹, α_(T), and α_(V), are determined before applying thegeomechanical equations to the modeling of a geomechanical reservoirsystem.

In certain examples, the reservoir permeabilities may be expressed as aninitial directional permeability (K_(abs)) multiplied by a permeabilitymultiplier ƒ for the permeability at every point and for every timestep: K=K_(abs)×ƒ, where ƒ is as function of one or more otherparameters, such as the fluid pressure, total stresses, bulk volumetricstrain, pore pressure, initial reference pressure, principal stresses,effective plastic strain, current porosity, initial porosity, anddeposition fraction.

In the constitutive equations for a transversely isotropic material,strains can be expressed in terms of stresses. In certain computations,effective stresses may be used for the calculations, and the initialin-situ stresses can be nonzero. In other computations, the equationsmay be simplified by using total stresses and using an assumption thatthe initial in-situ stresses can be zero. The constitutive equationsrelating stresses and strains for a transversely isotropic material canbe expressed as:

$\begin{matrix}{ɛ_{xx} = {{\frac{1}{E_{h}}\sigma_{xx}} - {\frac{v_{h}}{E_{h}}\sigma_{yy}} - {\frac{v_{v}}{E_{v}}\sigma_{zz}}}} & (17) \\{ɛ_{yy} = {{\frac{1}{E_{h}}\sigma_{yy}} - {\frac{v_{h}}{E_{h}}\sigma_{xx}} - {\frac{v_{v}}{E_{v}}\sigma_{zz}}}} & (18) \\{ɛ_{zz} = {{\frac{1}{E_{v}}\sigma_{zz}} - {\frac{v_{h}}{E_{v}}\sigma_{xx}} - {\frac{v_{v}}{E_{v}}\sigma_{yy}}}} & (19) \\{\gamma_{xy} = {\frac{2\left( {1 + v_{h}} \right)}{E_{h}}\sigma_{xy}}} & (20) \\{\gamma_{xz} = {\frac{1}{G_{v}}\sigma_{xz}}} & (21) \\{\gamma_{yz} = {\frac{1}{G_{v}}\sigma_{yz}}} & (22)\end{matrix}$

Equations 17-22 assume that the axis of symmetry is the z-direction andthat the vertical direction is the z-direction. Five elastic constantsin Equations 17-22 can be supplied to the model before a computationinvolving transversely isotropic solids is performed. It can be assumedthat Biot's constants are isotropic and that the thermal expansioncoefficients are isotropic when performing poroelastic computations, sotwo Biot constants and two thermal expansion coefficients may besupplied in addition to the five transversely isotropic constants whenanalyzing a transversely isotropic porous material. For a poroelasticcomputation including the thermal model with transversely isotropicelastic constants, the stresses in Equations 17-22 can be effectivestresses and the strains in Equations 17-22 can be effective strainsdefined as:

$\begin{matrix}{ɛ_{ij}^{e} = {ɛ_{ij} - {\frac{1}{3}{\alpha_{T}\left( {T - T_{o}} \right)}\delta_{ij}}}} & (23)\end{matrix}$

where ε_(ij) is the true strain.

5.5.2.2 Poroplastic Materials

A poroplastic material exhibits nonlinear behavior, in that it mayundergo permanent (i.e., plastic) volumetric strains, and hence,porosity changes. Large fluid pressures can cause the poroplasticmaterial to yield. As a result, geomechanical computations for aporoplastic material may predict large sudden changes in fluidporosities. These large sudden changes in porosity can cause significantstability problems and also produce negative fluid pressures. Thenegative pressures normally arise when the fluid compressibility is low,the permeability is low, and the porosity expansion is sudden and large.The equations for a poroelastic material discussed in Section 5.5.2.1above are applicable to poroplastic materials also. However, theporosity can be modified to account for the changes in porosity that maybe predicted for poroplastic materials.

In certain examples, an equation can be used to damp sudden porositychanges in order to improve the numerical stability of the computationsand to reduce the frequency of encountering negative pressures. In thesecomputations, the porosity in the reservoir model can be defined as afluid porosity (φ_(fluid)), and can be treated as different from theporosity in the geomechanical model (φ_(geomech)). The relationshipbetween fluid porosities and geomechanical porosities would be governedby the following damping equation during the computations:

$\begin{matrix}{{\tau \frac{\partial\phi_{fluid}}{\partial t}} = {\phi_{geomech} - \phi_{fluid}}} & (24)\end{matrix}$

The fluid porosity (φ_(fluid)) can be computed using Equation 24 andused in the fluid flow equations, while the geomechanical porosity(φ_(geomech)) can be computed in the geomechanical model, and τ is aprescribed time constant. After a step change in (φ_(geomech)), therelative difference between φ_(fluid) and φ_(geomech) can be less than1% after a time period of 5τ. The value of τ can be chosen so that thecomputations are stable and can be chosen to be as short a time intervalas possible, for example, by setting τ to a value that is shorter thanthe time frame of a time step in the computation or the total time ofthe computation. By way of example, if a computation is expected to lastseveral days, then τ can be on the order of minute; if a computation isexpected to last a few minutes, then τ can be on the order ofmilliseconds. The value of τ can be set to about one minute forcomputations on many poroplastic materials. In other examples, the valueof τ can be set to zero (which removes all damping). The value of τ touse for a given computation will be apparent to one skilled in the art.

5.5.2.3 Plastic Flow

Plastic flow, also called yielding, denotes a permanent deformation of amaterial. Once yielding is initiated in a material, plastic flow may ormay not persist. A yield condition is a mathematical representationwhich marks the transition from elastic to plastic deformation. Anassumption in the plastic flow equations is that a single yieldcondition and a single plastic potential can be used to describe thereservoir material, and that a single hardening parameter can be used torepresent the movement of the yield surface. It can be further assumedthat the stress and strain rates can be expressed as:

{dot over (σ)}_(ij) =E _(ijkl)({dot over (ε)}_(kl)−{dot over (ε)}_(kl)^(p))  (25)

where tensor E_(ijkl) denotes the linear elastic properties of thereservoir material:

E _(ijkl)=λδ_(ij)δ_(kl)+μ(δ_(ik)δ_(jl)+δ_(il)δ_(jk))  (26)

where λ and μ are Lame constants. Eq. 26 can be used for an isotropicelastic material. In certain computations, the plastic model may berestricted to isotropic elastic properties. In other computations, Eqs.25 and 26 may include other terms which model plastic properties of amaterial.

An assumption can be made that there exists a plastic potential G whichmay be equal to the yield criterion in some cases, but which may differfrom the yield condition for non-associated flow. Further, it can beassumed that the plastic strain rates are given by:

$\begin{matrix}{{\overset{.}{ɛ}}_{ij}^{p} = {\overset{.}{\lambda}\frac{\partial G}{\partial\sigma_{ij}}}} & (27)\end{matrix}$

where the value of the scalar {dot over (λ)}≧0, which is not related tothe Lame constant, is related to certain constraints which may be placedon the properties of the material, i.e., the yield condition. Theplastic potential G can be used to determine the directions of theplastic strain rates while the scalar can be used to determine theirmagnitudes.

An effective plastic strain (ε^(p)) can be defined as:

ε_(p)=∫{dot over (ε)}^(p) dt=∫√{square root over ({dot over (ε)}_(ij)^(p){dot over (ε)}_(ij) ^(p))}dt  (28)

and the relationship between the rate of change of the effective plasticstrain and the {dot over (λ)} parameter is:

$\begin{matrix}{{\overset{.}{ɛ}}^{p} = {\overset{.}{\lambda}\sqrt{\frac{\partial G}{\partial\sigma_{ij}}\frac{\partial G}{\partial\sigma_{ij}}}}} & (29)\end{matrix}$

If a unit tensor n_(ij) is defined as:

$\begin{matrix}{n_{ij} = \frac{\frac{\partial G}{\partial\sigma_{ij}}}{\sqrt{\frac{\partial G}{\partial\sigma_{kl}}\frac{\partial G}{\partial\sigma_{kl}}}}} & (30)\end{matrix}$

the plastic strain rate tensor may be written as

{dot over (ε)}_(ij) ^(p)={dot over (ε)}^(p) n _(ij)  (31)

It can be assumed that a simplified form of the yield condition can bewritten in terms of the stresses (σ_(ij)) and a hardening parameter (κ)as:

F(σ_(ij),κ)≦0  (32)

where F is negative in the elastic state and zero in the plastic state,and κ is a function of the plastic strains (strain-hardening). Duringplastic deformation, the yield condition satisfies the relation:

$\begin{matrix}{{{\frac{\partial F}{\partial\sigma_{ij}}{\overset{.}{\sigma}}_{ij}} + {\frac{\partial F}{\partial\kappa}\frac{\partial\kappa}{\partial ɛ_{ij}^{p}}{\overset{.}{ɛ}}_{ij}^{p}}} = 0} & (33)\end{matrix}$

Equation 33 may be combined with Equations 25, 26 and 27 to arrive atthe expression:

$\begin{matrix}{{\overset{.}{\sigma}}_{ij} = {\left\{ {E_{ijkl} - \frac{E_{ijpq}\frac{\partial G}{\partial\sigma_{pq}}\frac{\partial F}{\partial\sigma_{mn}}E_{mnkl}}{{\frac{\partial F}{\partial\sigma_{rs}}E_{rsab}\frac{\partial G}{\partial\sigma_{ab}}} - {\frac{\partial F}{\partial\kappa}\frac{\partial\kappa}{\partial ɛ_{ab}^{p}}\frac{\partial G}{\partial\sigma_{ab}}}}} \right\} {\overset{.}{ɛ}}_{kl}}} & (34)\end{matrix}$

The constitutive equation relates the stress rates to the strain rateswhere the coefficients depend on the elastic properties, the currentstresses, the current plastic strains, and a hardening parameter. For anassociative-type flow computation, it can be assumed that F=G, whichmakes the constitutive equation symmetric. If the nonlinear flowequations are solved using an implicitly-coupled computation, theJacobian for the system of equations may be slightly non-symmetric, evenwhen the constitutive equations are symmetric.

5.5.2.4 Hardening Models

An examination of the requirements for subsequent plastic deformationand the stress-strain relationship of the material can provide anindication of whether the plastic flow may or may not persist. Amaterial can be ideally plastic or subject to strain-hardening. Anideally plastic material (such as but not limited to structural steel)exhibits a yield condition which remains unaltered by plasticdeformation. However, many materials are altered by inelasticdeformation (termed strain-hardening) and the yield condition can bemodified as the materials' resistance to yielding increases.

The geomechanical model includes hardening models which can be used tomodel the plastic deformation of material within the geomechanicalreservoir system. Examples of hardening models include theDrucker-Prager model with shear hardening or cap hardening, the modifiedDrucker-Prager model with tabular hardening, the modifiedBigoni-Piccolroaz model, and the Matsuoka-Nakai model.

Drucker-Prager with Shear Hardening

The yield condition for the Drucker-Prager equation with shear hardeningand positive tensile stress can be expressed as:

F=√{square root over (J₂)}+αI ₁ −k(ε^(p))^(m)−Γ=0  (35)

where all scalar parameters are nonnegative constants, α is a constant,m is an exponent, and Γ is the yield constant. Values for the parametersα, k, m, and Γ correspond to the input parameters may be entered toperform a computation. The constant α may take on values between 0.0 and1/√3. If the value of α is zero, then the Drucker-Prager model becomes aVon Mises model of plasticity. The parameter α is related to thefriction angle that is used for a Mohr-Coulomb model. The effectiveplastic strain ε^(p) is a hardening parameter which may be computed andreported, for example, to a user.

After hardening, i.e., an increase in ε^(p), the yield surface may beconsidered to move from an original surface position to a final surfaceposition, as shown in FIG. 11A. In this computation, the first invariantof stress I₁ (the first invariant of total stress) is negative incompression (I₁=σ_(kk), summation over k).

The plastic flow equation for the Drucker-Prager model can berepresented:

$\begin{matrix}{{\overset{.}{ɛ}}_{ij}^{p} = {\overset{.}{\lambda}\left\{ {{\alpha \; \delta_{ij}} + \frac{S_{ij}}{2\sqrt{J_{2}}}} \right\}}} & (36)\end{matrix}$

where {dot over (ε)}^(p)={dot over (λ)}√{square root over (3α²+½)} and αin Equation 35 is the same α used in the yield equation if associatedflow is assumed.

Drucker-Prager with Cap Hardening

A Drucker-Prager models with cap hardening can have two yield surfaces.One yield surface is the non-hardening Drucker-Prager failure surface(perfectly plastic) given by:

F _(s) =√{square root over (J₂)}+αI ₁−Γ=0  (37)

where all scalar parameters are nonnegative constants and tensilestresses are positive, α is a parameter and Γ is the yield constant. Thesecond yield surface is an elliptical hardening cap of the form:

$\begin{matrix}{F_{c} = {{J_{2} - {\frac{1}{R^{2}}\left\lbrack {\left( {X - L} \right)^{2} - \left( {I_{1} - L} \right)^{2}} \right\rbrack}} = 0}} & (38)\end{matrix}$

The variable X≦RΓ is a function of the plastic volumetric strain(compaction is negative) and is the value of the first invariant ofstress I₁ at which the elliptical cap intersects the axis of I₁. Thevariable L≦0 is also a function of the plastic volumetric strain and isthe value of I₁ at which the elliptical cap intersects theDrucker-Prager surface F_(s). The elliptical cap is vertical at theintersection with the I₁ axis and horizontal at the intersection withthe Drucker-Prager failure surface. The following equations can be usedto relate X to the effective plastic volumetric strain and theconstraint enforcing zero slope for the cap at the Drucker-Prager yieldsurface:

$\begin{matrix}{{X = {X_{o} + {\frac{1}{D}{\ln \left( {1 + \frac{{\overset{\_}{ɛ}}_{v}^{p}}{W}} \right)}}}}{or}{X = {X_{o} - {H\left( {- {\overset{\_}{ɛ}}_{v}^{p}} \right)}}}{and}} & (39) \\{\frac{X - L}{R} = {{\alpha \; L} - \Gamma}} & (40)\end{matrix}$

Two cap hardening models can be derived from Equations 39 and 40. Thefirst cap hardening model uses the logarithm in Equation 39 forhardening at the cap. The second cap hardening model uses the tabularfunction H( ) of Equation 39, where H( ) is strictly monotonicallyincreasing and H(0)=0. The effective plastic volumetric strain, ε _(v)^(p), and the variables L and X are all zero or negative in Equations 39and 40. The initial value of ε _(v) ^(p) can be nonzero for acomputation when the initial magnitude of X exceeds the magnitude ofX_(o).

Movement of the yield surface is shown in FIG. 11B, where the solid lineis the original surface and the dashed line is the location of thesurface after hardening, and L₁ and L₂ are the values of L for theinitial and subsequent yield surfaces.

Using the plastic equations, it can be shown that the plastic volumetricstrain rate can be expressed in terms of the {dot over (λ)} elasticparameter from the equation {dot over (ε)}_(ij) ^(p)={dot over(λ)}∂F_(c)/∂σ_(ij) in the form

$\begin{matrix}{{\overset{.}{ɛ}}_{v}^{p} = {\frac{6}{R^{2}}\left( {I_{1} - L} \right)\overset{.}{\lambda}}} & (41)\end{matrix}$

where I₁≦L≦0 whenever plastic deformation occurs on the cap, and {dotover (ε)}_(v) ^(p)≦0 and {dot over (λ)}≧0. Hardening parameters whichmay be computed and reported, for example, to a user, are L, {dot over(ε)}_(v) ^(p), and {dot over (λ)}. The plastic flow equation on theshear surface can be the same as described previously for the standardDrucker-Prager model, while flow on the cap can be given by:

$\begin{matrix}{{\overset{.}{ɛ}}_{ij}^{p} = {\overset{.}{\lambda}\left\{ {{\frac{2}{R^{2}}\left( {I_{1} - L} \right)\delta_{ij}} + S_{ij}} \right\}}} & (42)\end{matrix}$

Modified Drucker-Prager with Tabular Hardening

The yield condition for the modified Drucker-Prager equation with shearhardening and positive tensile stress is

$\begin{matrix}{F = {{{\frac{1}{2}{\sqrt{J_{2}}\left\lbrack {1 + \frac{1}{K} - {\left( {1 - \frac{1}{K}} \right){\cos \left( {3\; \theta} \right)}}} \right\rbrack}} + {\alpha \; I_{1}} - {H\left( ɛ^{p} \right)} - \Gamma} = 0}} & (43)\end{matrix}$

where α is a constant, K is the deviation, Γ is the yield constant, andθ is π/3 for a uniaxial compression test and can be defined by

$\begin{matrix}{\theta = {\frac{1}{3}{\cos^{- 1}\left( {\frac{3\sqrt{3}}{2}\frac{J_{3}}{J_{2}^{3/2}}} \right)}}} & (44)\end{matrix}$

H( ) is a tabular function of the effective plastic strain, I₁ is thefirst stress invariant, and J₂ and J₃ are the second and thirdinvariants of the total stress deviator:

$\begin{matrix}{J_{2} = {\frac{1}{2}S_{ij}S_{ij}}} & (45) \\{J_{3} = {\frac{1}{3}S_{ij}S_{jk}S_{ki}}} & (46) \\{S_{ij} = {\sigma_{ij} - {\frac{1}{3}\sigma_{kk}\delta_{ij}}}} & (47)\end{matrix}$

where S_(ij) is the deviatoric stress (minus the mean stress), δ_(ij) isthe kronecker delta (which is 1 when i=j and 0 otherwise), J₂ is summedover i and j, and J₃ is summed over i, j, and k. The modifiedDrucker-Prager model reduces to the Drucker-Prager model if the materialparameter K is set to 1.0 and if the table values of H( ) satisfyH(ε^(p))=k(ε^(p))^(m). The value of parameter K may range from 0.78 to1.0. The hardening parameter of effective plastic strain ε^(p) may bereported, for example, to a user. The constants α, K, and Γ may beentered.

The effect of the parameter K is illustrated in FIGS. 12A-D, which showplots of the modified Drucker-Prager yield surface in the octahedral(deviatoric) plane for four values of K in the octahedral plane (and fora constant I₁). In the four plots of FIGS. 12E-H, it can be assumed thatcompression is positive, and all have been normalized to cross 1.0 onthe positive vertical axis. The standard Drucker-Prager has K=1.0 and iscircular in the plane. The modified Drucker-Prager surface is no longerconvex for values of K less than 0.78. Certain computations may berestricted so that values of K are greater than or equal to 0.78.

The plastic flow equation for the Modified Drucker-Prager equation is

$\begin{matrix}{{\overset{.}{ɛ}}_{ij}^{p} = {\overset{.}{\lambda}\left\{ {{\alpha_{f}\delta_{ij}} + {\frac{S_{ij}}{4\sqrt{J_{2}}}\left\lbrack {1 + \frac{1}{K} - {\left( {1 - \frac{1}{K}} \right){\cos \left( {3\; \theta} \right)}}} \right\rbrack} + {\frac{3\sqrt{3}}{4}\left( {1 - \frac{1}{K}} \right)T_{ij}}} \right\}}} & (48)\end{matrix}$

where T_(ij) is defined by

$\begin{matrix}{T_{ij} = {{\frac{2}{3}\delta_{ij}} - \frac{S_{ik}S_{kj}}{J_{2}} + {\frac{\cos \left( {3\; \theta} \right)}{\sqrt{3\; J_{2}}}S_{ij}}}} & (49)\end{matrix}$

and α_(f) in the flow equation is the same used in the yield equationassociated flow is assumed. If it is assumed that T_(ii)=0 and S_(ii)=0,setting α_(f) equal to zero in Equation 49 results in the plasticvolumetric strains being zero for a computation.

Modified Bigoni-Piccolroaz with Tabular Hardening

The modified Bigoni-Piccolroaz model (MBP) is a modification of theyield criteria in D. Bigoni et al. (2004) “Yield criteria forquasi-brittle frictional materials,” Intl. J. of Structures41:2855-2878. The yield equation for the MBP model with positive tensilestress can be expressed as:

$\begin{matrix}{F = {{{A\; \cos \; (\phi)\sqrt{J_{2}}} + {\alpha \; I_{1}} - {H\left( ɛ^{p} \right)} - \Gamma} = 0}} & (50) \\{A^{- 1} = {\cos \left\{ {{\frac{1}{3}{\cos^{- 1}\left( {- \gamma} \right)}} - {\beta \frac{\pi}{6}}} \right\}}} & (51) \\{\phi = {{\frac{1}{3}\cos^{- 1}\left\{ {\gamma \; {\cos \left( {3\; \theta} \right)}} \right\}} - {\beta \frac{\pi}{6}}}} & (52)\end{matrix}$

where α is a constant, Γ is the yield constant, γ is the deviation, andwhere 0≦α≦1/√{square root over (3)}, 0≦β≦1, Γ≧0, and 0≦γ1 are materialconstants and H( ) is a tabular function of the effective plasticstrain. Parameters J₂ and J₃ are invariants of the stress deviator, andI₁ is the first stress invariant. The constant A can be chosen such thatA cos(φ) is 1.0 for a uniaxial compression test. The MBP model usesstrain-hardening to model plastic hardening and an assumption can bemade that the yield surface expands uniformly (no rotation) as ithardens. The constants α, β, Γ, and γ, which can be input parameters,control the shape of the yield surface in principal stress space. Theangle of the yield cone in principal stress space is controlled by α,while the location of the tip of the cone on the I_(i) axis is given byΓ/α before hardening. The shape of the yield surface in the octahedralplane (for constant I₁) can be controlled by the parameters β and γ.Several common yield surfaces can be reproduced for specific choices ofthe MBP parameters. Common yield surfaces that are special cases of theMBP yield surface are von Mises, Drucker-Prager, Mohr-Coulomb, Lade,Tresca, and Matsuoka-Nakai.

Parameters β and γ can be modified to determine the shape of the yieldsurface in the octahedral plane. The yield surfaces which result fordifferent values of β and γ are shown in FIGS. 12E-H. In the plots ofFIGS. 12E-H, it can be assumed that compression is positive, and allplots have been normalized to cross 1.0 on the positive vertical axis.The MBP model yield equation does not exhibit the loss of convexitywhich can occur for the modified Drucker-Prager model yield equation.For example, for the case where β=0.0 and γ=1.0, the MBP model yieldequation generates a triangle in the octahedral plane (as compared tothe non-convex surfaces for the modified Drucker-Prager model yieldequation when K=0.60).

The plastic flow equation for the MBP model is

$\begin{matrix}{{\overset{.}{ɛ}}_{ij}^{p} = {\overset{.}{\lambda}\left\{ {{\alpha_{f}\delta_{ij}} + \frac{A\; \cos \; (\phi)S_{ij}}{2\sqrt{J_{2}}} - {\frac{A\; \gamma \sqrt{3}{\sin (\phi)}}{2\; {\sin \left( {{3\; \phi} + {\beta \frac{\pi}{2}}} \right)}}T_{ij}}} \right\}}} & (52)\end{matrix}$

and α_(f) in the flow equation is the same α used in the yield equationassuming associated flow. If α_(f) is set equal to zero in Equation 52,the plastic volumetric strains is zero for the computations.

5.5.2.5 Sand Production Model

One or more parameters from a hardening model for computations of theplastic deformation may be used for computation of the sand productionmodel. In computations using the sand production model, at least onesand production criterion can be applied to one or more grid cells. Itmay be assumed that a computation grid cell fails when a critical valueis reached for:

(1) the total strain invariant, or

(2) the plastic strain invariant, or

(3) the maximum effective stress

at the center of a grid cell. With each of these sand productioncriteria, failure of a grid cell represents sanding of reservoirmaterial (resulting in a cavity).

For criterion (1), the total strain invariant can be expressed as:

ε=√{square root over (ε_(ij)ε_(ij))}  (53)

For criterion (2), the plastic strain invariant can be expressed as:

ε^(p)=√{square root over (ε_(ij) ^(p)ε_(ij) ^(p))}  (54)

In certain computations, the plastic strain invariant criterion may notbe applied if the hardening model is the Drucker-Prager cap hardeningmodel.

For criterion (3), the maximum effective stress, the maximum principaleffective stress can be computed and its value can be compared to theinput value (where tensile stresses are assumed positive).

Criterion (2), the plastic strain invariant can be used to account fortransient sand production from a cell prior to the total failure of thatcell. If the critical plastic strain for total cell failure is ε_(lim)^(p), then the fraction of sand produced from a cell prior to total cellcollapse can be represented by a production function. The productionfunction can be a function ƒ(x), where ƒ(x)=0 when x=0, and ƒ(x)=1 whenx=1, and where x is a function of the critical plastic strain(x=g(ε_(lim) ^(p))). The function ƒ(x) may be any monotonic function ofx, including but not limited to, a fraction function, a power function,a sine function, a cosine function, a logarithmic function, anexponential function, a sigmoid function, or any combination thereof. Insome examples, x can be a function of both the critical plastic strain(ε_(lim) ^(p)) and the plastic strain invariant (ε^(p)), such as but notlimited to a ratio ε^(p)/ε_(lim) ^(p). In an example, the productionfunction can be given by:

$\begin{matrix}\left( \frac{ɛ^{p}}{ɛ_{\lim}^{p}} \right)^{m} & (55)\end{matrix}$

where m is an exponent.

5.5.2.6 Computation in the Geomechanical Model

Computation of the geomechanical model can be on the 3D grid (which canbe the grid used for the reservoir model). For example, a standardfinite element method may be used for the geomechanical equations wherestresses are integrated at the eight Gaussian points interior to eachelement, and fluid pressures are integrated at the center of eachelement. The discretization produces a 27-point stencil for thedisplacements when eight elements share a common corner, and there arethree displacements at each node. In certain computations, thedisplacements are the primary variables for the geomechanical equationswhere the displacements are evaluated at the corners (node-based) of thehexahedral elements.

5.5.3 Thermal Model

The thermal model calculates temperature changes that occur when hot orcold fluids are injected into or produced from a reservoir and alsocalculates conduction from a well-bore that may be circulating hot/coldfluids but is not actually injecting or producing fluids. Thetemperature calculations may include porous flow and geomechanics butmay be configured to not include steam injection. Injected water can beassumed to be in a liquid phase.

The thermal model calculates temperature changes that arise due toconduction and convection in the reservoir, hot or cold fluidinjection/production at wells, conduction at wells, and thermal andmechanical interactions between the fluids and solid. These mechanismscan be combined in a single energy equation that is solved along withthe porous flow equations and solid deformation equations. The energyequation can be formulated in terms of a Lagrangian-approach for theporous solid and all fluid movement can be relative to the movement ofthe porous solid.

With this type of formulation, the mass of the porous solid can beconstant when evaluating the energy balance for an element of thereservoir while the mass of fluid changes as fluids flow into and out ofthe porous solid.

5.5.3.1 Combined Accumulation Term

The combined expressions for energy change in the fluid and solid (whenthe geomechanical system is included in the computation) can be givenby:

$\begin{matrix}{{\overset{.}{u}}_{T} = {{\frac{\partial\;}{\partial t}\left( {\phi {\sum\limits_{\alpha}^{\;}\; {\rho_{\alpha}S_{\alpha}h_{\alpha}}}} \right)} + {C_{r}\overset{.}{T}} + {\alpha_{T}{KT}_{o}{\overset{.}{ɛ}}_{ii}} - {\left( {{\alpha_{V}T_{o}} + \phi_{o}} \right)\overset{.}{p}} + {\sigma_{ij}{\overset{.}{ɛ}}_{ij}}}} & (56)\end{matrix}$

where φ is the porosity with respect to the undeformed bulk volume, α isthe phase, ρ_(α) is the density of phase, S_(α) is the saturation ofphase α, h_(α) is the specific enthalpy of phase α, α_(T) is thevolumetric expansion coefficient in stress/strain equations, K is thepermeability, α_(v) is volumetric expansion coefficient in the porosityequation, and T₀ is the fluid temperature. Conduction and convectionterms can affect the combined accumulation term. In Equation 56, theapproximation can be made that the phase pressures are the same, theterm φ{dot over (p)} can be approximated by φ_(o){dot over (p)}, and theterm C_(r) can be C_(r)≡C_(v)+α_(V) ²MT_(o). The heat capacity C_(v) isthe heat capacity for a porous solid measured at constant pore volumeand constant bulk volume. It can be shown that the heat capacity C_(r)as defined above is the heat capacity at constant fluid pressure andconstant bulk volume. Equation 56 is the general formula for the totalrate of change of internal energy in an element and is used for theaccumulation term when terms including the fluid temperature T₀ areconsidered in the computation and the computation includes both thermaland geomechanical models.

Equation 56 takes on simpler forms if the fluid temperature T₀ isconsidered in the computation.

For computations that include thermal and geomechanical models and termsincluding the fluid temperature T₀ and φ_(o){dot over (p)} are notconsidered in the computation, the simplified expression is

$\begin{matrix}{{\overset{.}{u}}_{T} = {{\frac{\partial\;}{\partial t}\left( {\phi {\sum\limits_{\alpha}^{\;}\; {\rho_{\alpha}S_{\alpha}h_{\alpha}}}} \right)} + {C_{r}\overset{.}{T}} + {\sigma_{ij}{\overset{.}{ɛ}}_{ij}}}} & (57)\end{matrix}$

For computations that include the thermal model, but do not include thegeomechanical model, the accumulation term is:

$\begin{matrix}{{\overset{.}{u}}_{T} = {{\frac{\partial\;}{\partial t}\left\lbrack {\phi\left( {{\sum\limits_{\alpha}^{\;}\; {\rho_{\alpha}S_{\alpha}h_{\alpha}}} - p} \right)} \right\rbrack} + {C_{r}\overset{.}{T}}}} & (58)\end{matrix}$

while for computations that include the thermal model, but do notinclude the geomechanical model, and the fluid internal energy isapproximated by the fluid enthalpy, the expression is:

$\begin{matrix}{{\overset{.}{u}}_{T} = {{\frac{\partial\;}{\partial t}\left( {\phi {\sum\limits_{\alpha}^{\;}\; {\rho_{\alpha}S_{\alpha}h_{\alpha}}}} \right)} + {C_{r}\overset{.}{T}}}} & (59)\end{matrix}$

5.5.3.2 Energy Equation

The choice of the accumulation term above determines the equations thatare included in the final energy conservation equation. The energyequation can be expressed as

$\begin{matrix}{{{\frac{\partial\;}{\partial t}\left( {\phi {\sum\limits_{\alpha}^{\;}\; {\rho_{\alpha}S_{\alpha}h_{\alpha}}}} \right)} + {C_{r}\overset{.}{T}}} = {{\nabla{\cdot \left( {K_{T}{\nabla T}} \right)}} - {\nabla{\cdot \left( {\sum\limits_{\alpha}^{\;}\; {\rho_{\alpha}h_{\alpha}{\overset{\rightarrow}{v}}_{\alpha}}} \right)}} + {\sum\limits_{\alpha}^{\;}\; {\rho_{\alpha}h_{\alpha}q_{\alpha}{\delta \left( \overset{\rightarrow}{x} \right)}}}}} & (60)\end{matrix}$

where a constant temperature boundary condition can be applied atinjection wells. Since Equation 60 does not contain any mechanicalterms, other than transport, the temperatures in grid blocks should notbe affected by the expansion or compression of a fluid phase or solidwhen the enthalpies are only functions of temperature. However, Equation60 does still account for heat of vaporization from a liquid to vaporbecause the latent heat can be included in the difference between theenthalpy of a component in a liquid phase and its corresponding enthalpyin the vapor phase.

If a simulation includes geomechanical calculations and terms includingthe fluid temperature T₀ are considered in the computation, then thegeneral energy equation can be given by:

$\begin{matrix}{{{\frac{\partial\;}{\partial t}\left( {\phi {\sum\limits_{\alpha}^{\;}\; {\rho_{\alpha}S_{\alpha}h_{\alpha}}}} \right)} + {C_{r}\overset{.}{T}} + {\alpha_{T}{KT}_{o}{\overset{.}{ɛ}}_{ii}} - {\left( {{\alpha_{V}T_{o}} + \phi_{o}} \right)\overset{.}{p}}} = {{\nabla{\cdot \left( {K_{T}{\nabla T}} \right)}} - {\nabla{\cdot \left( {\sum\limits_{\alpha}^{\;}\; {\rho_{\alpha}h_{\alpha}{\overset{\rightarrow}{v}}_{\alpha}}} \right)}} + {\sum\limits_{\alpha}^{\;}\; {\rho_{\alpha}h_{\alpha}q_{\alpha}{\delta \left( \overset{\rightarrow}{x} \right)}}}}} & (61)\end{matrix}$

The σ_(ij){dot over (ε)}_(ij) term is present on both sides of theenergy equation; i.e., in the accumulation term and also as work done onthe bulk solid, and consequently σ_(ij){dot over (ε)}_(ij) is notincluded in Equation 61.

5.5.3.3 Modifications to Fluid Densities and Viscosities

There are several options for computing temperature-dependent densitiesand viscosities for fluids, and the methods available for computingdensities and viscosities vary depending on the phase behavior model andthe numerical technique that is being used. In certain examples, thefluid properties can be computed directly as functions of temperatureand pressure. In certain other examples, the fluid properties can becomputed as a function of the current cell pressures and initial celltemperatures (isothermal flash) and modification functions can be usedto correct for the differences between the initial cell temperatures andthe current cell temperatures.

Fluid properties may be computed directly as a function of temperatureand pressure for implicit single phase runs, implicit two-phase runs,and implicit compositional runs. These computations involving thegeomechanical model can be iteratively coupled or fully coupled.

The viscosities and densities that are calculated in the black-oil,K-value, and compositional PVT packages may be modified to account forchanges in fluid properties due to temperature changes. The initialreservoir temperature and current fluid pressure during a computationmay be used to calculate the fluid properties in each of the PVTpackages (isothermal flash) and then these properties are modified bymultiplying by user-supplied modification factors, i.e.,μ(p,T)=μ(p,T_(o)) μ(T) where μ(T) can be entered as a table andμ(p,T_(o)) can be calculated in the PVT package.

For many thermal-geomechanical studies, interest can be primarilyfocused on how temperature affects fluid viscosities and thermalstresses. For these types of applications, either the modificationfactors or direct calculation of fluid properties may be used, and theresults can be very similar using both techniques. However, forapplications where temperature has a strong effect on densities and/orphase separation, then the direct calculation of fluid properties can beused as functions of temperature and pressure.

5.6 Interrelationship of the Various Models

A system and method can couple the interactions between the differentgeomechanical reservoir system models, e.g., between porous flow,geomechanical displacements and stresses, sand production and cavityformation. For example, the geomechanical solution can influence theporous flow computations through the porosity and permeability terms. Asanother example, the fluid solution in the reservoir can affect thegeomechanical calculations through its role in the effective stresses.In yet another example, the fluid solution can affect the geomechanicalcalculations involving the sand production model through normal tractionthat occurs at the face of the cavity and can affect the flow in thereservoir.

5.7 Implementation of the Sand Production Prediction

5.7.1 Computations Including the Sand Production Model

FIG. 13 illustrates several considerations in an example computation ofthe sand production prediction in a time step of a computation where acavity has been generated between the reservoir overburden (rockoverlying the area of interest in the well) and the underburden (rockunderlying the area of interest in the well). At the end of each timestep, a check can be made for failed grid cells, and a cell can beremoved immediately from the deformation computations (for example, acomputation involving a geomechanical model) once it is identified asbeing a failed cell. A failed cell is a cell for which a sand productioncriterion is reached (see Section 5.5.2.5 supra), which can be used toindicate the point where the material fails, i.e., rubblizes to formsand and generate a cavity. Once a cell fails and becomes part of thecavity, it can be excluded from the stiffness matrix for computing thedisplacements. As indicated in FIG. 13, unfailed reservoir cells aheadof the interface between an existing cavity and unfailed cells can bechecked for failure in each time step. As also illustrated in FIG. 13,the traction condition at the interface between an existing cavity andunfailed cells can be taken to be equal to the fluid pressure in thecavity cells at interface. The standard porous flow models andgeomechanical models can be applied to unfailed, intact grid cells.

The volume of sand in the failed cell can be then added to the reportedsand production volume for the nearest production well. In somecomputations, the sand volume of a cell can be the initial bulk volumeof that cell. The failed cell can be retained in the fluid flowcalculations to allow flow between the wall of the generated cavity andthe well, and the porosity of the cell can be kept constant. In somecomputations, the permeabilities can be increased in the cavity tominimize the pressure drop between the well and the wall of the cavity(for fluid flow). Once a cell is removed from the deformationcalculations, its properties no longer affect the deformations, exceptthat the failed cells immediately adjacent to the cavity can exert acompressive load on the unfailed porous solid (i.e., an unfailed cell)and this compressive load at the wall of the cavity can be directlyrelated to the fluid pressures in the failed cells adjacent to the wallof the cavity.

Cavity generation calculations can be physically and/or numericallyunstable. To avoid physical instabilities, some hardening can be addedto the end of a curve that is nearly elastic-perfectly plastic, and atime frame can be set for the growth of an unstable cavity (i.e., thegrowth rate for the cavity can be directly related to the time scaleset, such as, as a non-limiting example, one minute. In order tominimize numerical instabilities, several parameters of the models maybe set to certain values. For example, the time constant τ from Equation24 (which relates changes in fluid porosities to changes ingeomechanical porosities) can be set to a non-zero value to minimizeinstabilities, since a value of zero removes all damping. A higher valueof time constant τ may enhance numerical stability. As another example,the permeability assigned to grid cells that fail (indicating sanding)can be set to a value which minimizes numerical instability. As anon-limiting example, the permeability can be set to a value on theorder of about 1000 Darcies. Lower or higher values may be set for thepermeability to enhance numerical stability.

5.7.2 Solution of the System of Partial Differential Equations

The Jacobian for the full system of equations for the models included inthe computation can be fully expanded and all equations can be solvedsimultaneously in the linear solver. The various mechanisms can becoupled using either one-way, explicitly, iteratively, or fully coupledtechniques. The most stable coupling technique can be used for sandproduction predictions.

The equations discussed in the previous sections can be implemented in asingle program and the program computes an implicit solution for all theequations, i.e., a backward-Euler technique can be used to approximatethe temporal derivatives and all primary variables and coefficients inthe equations can be evaluated at the end of the time step. Non-limitingexamples of primary variables are the fluid pressures in the reservoir,the displacements in the reservoir and at the boundary of the reservoir(if unconstrained), the fluid pressures in a fracture (or fluid volumesin the fracture if partially saturated), the well pressure, and at leastone parameter of the production function. The system of equations can besolved using a Newton-Raphson technique where the fully-expandedJacobian is generated for the entire system of equations and incrementalcorrections are found using a linear solver that includes all thesolution variables.

Cavity generation can be simulated by computing flow to a failed celland computing the effective stresses and displacements that arise in theadjoining reservoir. Sand production criteria are then combined with thestress and displacement states and traction condition at the cavityinterface to decide if the cavity has progressed. The cavity mayprogress across multiple grid cells during a single timestep. Eachtimestep can begin with the cavity configuration, geomechanical state,fluid state (if applicable), and thermal state (if applicable) from theprevious timestep. The computation then iterates to a converged solutionfor the conservation of mass and stress equilibrium equations assumingthe cavity does not propagate. Once a new converged solution isobtained, the sand production criteria can be checked to see if thecritical value of the sand production criterion has been reached in anycells. These cells can be treated as described in Section 5.7.1. Thissequence of iterations may be continued on the equations and checkingsand production criteria for a number of times, or until no furtherprogression of the cavity is seen for a timestep before moving on to thenext timestep.

5.8 Examples of Apparatus and Computer-Program Implementations

The methods disclosed herein can be implemented using an apparatus,e.g., a computer system, such as the computer system described in thissection, according to the following programs and methods. Such acomputer system can also store and manipulate the data indicative ofphysical properties associated with a geomechanical reservoir system,the fully-expanded Jacobian for the full system of equations for themodels included in the computation, the solution to the fully-expandedJacobian, the generated fracturing predictions, or measurements that canbe used by a computer system implemented with the analytical methodsdescribed herein. The systems and methods may be implemented on varioustypes of computer architectures, such as for example on a single generalpurpose computer, or a parallel processing computer system, or aworkstation, or on a networked system (e.g., a client-serverconfiguration such as shown in FIG. 14).

As shown in FIG. 14, the modeling computer system to implement one ormore methods and systems disclosed herein can be linked to a networklink which can be, e.g., part of a local area network (“LAN”) to other,local computer systems and/or part of a wide area network (“WAN”), suchas the Internet, that is connected to other, remote computer systems.

The modeling system comprises any methods of the described herein. Forexample, a software component can include programs that cause one ormore processors to implement steps of accepting a plurality ofparameters indicative of physical properties associated with thegeomechanical reservoir system, and/or the generated fracturingpredictions and storing the parameters indicative of physical propertiesassociated with the geomechanical reservoir system, and/or the generatedfracturing predictions in the memory. For example, the system can acceptcommands for receiving parameters indicative of physical propertiesassociated with the geomechanical reservoir system, and/or parameters ofa generated fracturing prediction, that are manually entered by a user(e.g., by means of the user interface). The programs can cause thesystem to retrieve parameters indicative of physical propertiesassociated with the geomechanical reservoir system, and/or parameters ofa generated fracturing prediction, from a data store (e.g., a database).Such a data store can be stored on a mass storage (e.g., a hard drive)or other computer readable medium and loaded into the memory of thecomputer, or the data store can be accessed by the computer system bymeans of the network.

6. RESULTS 6.1 Sand Production Prediction Tests

Sanding tests were performed on brine and kerosene saturated coresamples. The sand production prediction was applied to the results ofthe sanding tests. Temperature effects were not considered. Anassumption was made that single phase simulations could be used tosimulate the reservoir system.

6.2 Data Collection

The geometry of the grid, information on material behavior, boundaryconditions and permeability function are discussed below.

6.2.1 Geometry/Grid

The samples tested were approximately 100 mm long and had a circulardiameter of about 100 mm. The central hole was roughly 20 mm indiameter. The exact dimensions of both the brine and the kerosene testsamples are given in Table 1. FIGS. 15A and 15B show the core samplesused for the brine test, including the dimensions of the sample. Duringthe test, end-platens were used to apply the axial load. A hollowcylinder test was applied (such as illustrated in FIG. 3).

An axisymmetric grid of cell was set up for the computations, based onthe dimensions given in Table 1. The rubber sleeve was not representedin the computation grid cells. However the simulation model was setup toinclude the two platens, whose the radial dimensions match those of thesample and whose axial dimension was 20 mm.

TABLE 1 Sample dimensions. Length Diameter Interior hole diameter Test(mm) (mm) (mm) Brine 102.72 99.77 19.56 Kerosene 105.25 99.95 19.47

6.2.2 Material behavior

Physical constants of the materials were determined using data from atriaxial test. The measured material response and that computed usingthe models are in close agreement for both the brine saturated cores(shown in FIG. 16) and the kerosene saturated cores (shown in FIG. 17).The physical constants used for modeling the brine saturated sample aregiven in Table 2.

TABLE 2 Physical properties for the brine saturated sample. VariableValue Porosity  0.296 Poisson's ratio  0.03 Young's modulus 218823 psiRock density   2.65 g/cc Yield constant Γ 43.2 Alpha  0.36581616

6.2.3 Boundary conditions

The simulations were performed by solving the fully coupledgeomechanical and fluid flow models, hence both fluid flow andgeomechanical boundary conditions needed to be set (shown in FIG. 18).For the geomechanical boundary conditions, constant stress was assignedto the top, outer and inner surface (boundary) of the hollow cylinder(shown by arrows in panel (A) of FIG. 18). The geomechanical model wasdisplacement constrained at the bottom to keep the numerical model frommoving and rotating (i.e., the entire surface is constrained not to movein the z direction). To simulate test fluid flow conditions (see panel(B) of FIG. 18), a pressure drop is applied over the hollow cylinder,where the pressure on the outside of the cylinder is greater than thepressure on the inside of the cylinder. No fluid flow boundaryconditions were applied at the top and the bottom of the cylinder. Theporosity of the spacers was set to zero, which effectively made thepermeability of the spacers equal to zero as well.

The pressure on the inside of the hollow cylinder as well as thestresses on the inside surface of the hollow cylinder were kept equal to0 psi. The simulation used absolute pressure values, therefore a zerovalue of pressure was set to be 14.7 psi to account for atmosphericpressure. Both the stresses on the outside and the pressures on theoutside of the cylinder were varied dependent on the flow and stressregime applied in the test performed on the brine saturated sample. FIG.19 shows the measured and discretized (simulated) confining stress, flowrate and flowing pressure vs. time for the brine saturated sample.

6.2.4 Permeability Function

Using the test measurements, the permeability for a radial isotropichomogeneous geometry was computed. Darcy's law for radial flow for thecore sample geometry can be expressed as:

$\begin{matrix}{k = \frac{Q\; \mu \; {\ln \left( \frac{r_{e}}{r_{i}} \right)}}{2\; \pi \; h\; \Delta \; P}} & (62)\end{matrix}$

where k is the permeability, Q is the flow rate, μ is the viscosity,r_(i) is the inner radius, r_(e) is the outer radius, h is the heightand ΔP is the pressure drop. The permeability was determined using thedimensions of the core samples and Equation 62. The permeabilitydistribution resulting from this calculation showed that thepermeability was not constant as a function of flow rate (as shown inFIG. 20). Specifically as the flow rate increased, the permeabilitydecreased (observed to be affected by the quality of the pressuremeasurements). An average permeability was determined by fitting apower-law curve to a plot of the computed permeability as a function ofconfining stress (FIG. 21). This power-law function was used to computethe pressure to achieve the flow rate. In other words, the pressureswere altered such that, given the power-law fit derived permeabilityfunction, the desired flow rate was obtained.

Generally, the permeability is a function of cell properties and not ofthe boundary conditions, such as the confining stress. The assumptionmade was that the mean stress can be used as a proxy for the confiningstress, where the mean stress was calculated using:

$\begin{matrix}{\sigma_{mean} = {\frac{1}{3}\left( {\sigma_{1} + \sigma_{2} + \sigma_{3}} \right)}} & (62)\end{matrix}$

where σ_(mean) is the mean stress and σ₁ through σ₃ are the principalstresses. This assumption caused the computed flow rate and the measuredflow rate not to be identical. The flow rate measured in the tests wascompared to the flow rate predicted by the numerical simulation. Acorrection factor was applied to the outer boundary pressures to makethe flow rates match. The calculated pressures were multiplied by 1.035to obtain the flow rate match (shown in FIG. 22), which meant that onlya 3.5% correction factor was used. This flow rate match was obtained fora sample without sand production. Flow rates varied depending on thevolume of sand produced, as discussed in Section 6.3 infra.

6.3 Analysis

The computed sand match and its sensitivity to certain parameters arediscussed below. Also discussed is a question of whether a sustainable,stable cavity can be created.

6.3.1 Sand Production Model Computation

The following inputs for the sand production model were computed:

-   -   1. value of a critical plastic strain—the critical value of        effective plastic strain that causes the rock to fail.    -   2. exponent—exponent of the production function used for        computing the amount of sand production from a cell prior to        reaching the critical plastic strain value.    -   3. parameters of hardening model—the triaxial test data was        matched (fit) using a version of the modified Bigoni-Piccolroaz        (MBP) hardening model.

To fit the MBP hardening model to the triaxial test data, besides usingthe parameters given in Table 2, the MBP model can be made to allow auser to define the shape of the yield surface in the octrahedral plane(constant I₁) (see FIGS. 23A and 23B). The shape of the yield surface inthe octahedral plane is defined by two parameters for the MBP hardeningmodel: deviation (γ) and beta (β). Both the deviation and beta areconstants having a value between 0 and 1.0. Setting beta to 0 and thedeviation to 0 in the MBP model produces a yield surface in theoctahedral plane of the Drucker-Prager material model (shown in FIG.23A). Setting beta to 0 and the deviation to 0.95 in the MBP modelproduces a yield surface in the octahedral plane of the Lade-Typematerial model (shown in FIG. 23B).

Yielding occurs when a sample reaches the yield surface. Hence, linearelastic behavior governs at the interior areas of the shapes shown inFIGS. 22A and 22B. From FIGS. 22A and 22B, it was seen that thelikelihood of a stress path hitting the yield surface in the deviatoricplane for the Lade-Type model is higher than that for the Drucker-Pragertype model. It was expected that more yielding would occur with higherdeviation values when applying the two material models to general stressstates. It was noted that the yield surfaces in FIGS. 23A and 23Bproduced the same results when applied to triaxial compression tests.

6.3.1.1 Definition of Grid Cells

Computations were performed on 2D axisymmetric grids. The grid cell sizefor representing the sample was varied to investigate grid sizedependencies on the model results. The cell size of the sample was 1.284mm in the vertical direction and was varied from 0.05 mm to 0.01 mm to0.005 mm in the horizontal direction. Measurements were accurate to theorder of 0.1 cc, hence small cell sizes were used in the model tocapture this level of precision. The cell size in the spacers was 2 mmin the vertical direction. This resulted in models with 80300, 401100and 802100 cells respectively. Runtimes were observed to progressivelyincrease from around 2 hours to ˜14 hours to approximately 30 hours.FIG. 24 shows the progression of the model results as a function of cellsize. No noticeable change was observed in the predicted sand productionvolume from the grid with a horizontal cell size of 0.01 mm to 0.05 mm.The reduce runtimes below 14 hours, an additional grid cell model wascreated which had the same vertical cell dimensions as given before,however only the first 7 mm of the sample was gridded with a cell sizeof 0.01 mm. Beyond the initial 7 mm, the cell size was increased to0.11035 mm. This graduated grid cell model had 100,000 cells and ran inapproximately 3˜4 hours. The computation with this grid cell model wascompared with the results from other grid cell models, and it was seenthat the results for this graduated grid cells model overlapped with thefiner grid cell models, indicating that the results are independent ofthe cell size (FIG. 24). This graduated grid configuration was used forall subsequent runs.

6.3.1.2 Fit to the Production Function

The production function used for the fit was ƒ=(ε^(p)/ε_(lim) ^(p))^(m),and the simulations were performed by solving the system of partialdifferential equations of the coupled geomechanical model (including theproduction function) and the fluid flow model. The best fit was obtainedfor a deviation of 0.5, a beta value of 0, an exponent (m) value of 2,and a critical strain limit of 0.017 (see FIG. 25). The data was fit forboth the massive sanding and the gradual sanding portions of the curve.The onset of plastic deformation is indicated in FIG. 25. The water ratematch is displayed in FIG. 26. Close matches of both the sanding volumeand water rate were observed. A slight mismatch in the water rate isobserved at later times, due to the large permeability values assignedto cells that have sanded. The test results indicated that thepermeability increase slightly overestimated the water flow rate in thenumerical model. Furthermore, it was observed that the critical strainvalue of 0.017 indicated that sanding occurred far beyond the initialyielding of the sample (FIG. 27), indicating that, when the initialyield point is taken as the onset of sanding, the sanding tendencies ofa sample might be significantly overestimated.

To obtain the fit to the produced sand volume, first the deviation wasvaried, where the exponent (m) was kept constant at a value of 1. Foreach value of deviation, the critical strain value was chosen toresemble the yield value at massive sanding. To obtain this value, thefit was performed using a geomechanical model which contained thehardening model but not the sand production model. Multiple values ofdeviation, critical strain and exponent were used to investigate theuniqueness of the solution.

First the deviation value was varied from 0.92 to 0.5. As stated before,the critical plastic strain value was determined for each deviation byperforming the fit using a geomechanical model which contained thehardening model but not the sand production model. The yield value atthe onset of massive sanding was used as the initial value of thecritical plastic strain. Adjustments to the value of the criticalplastic strain were sometimes made to obtain the best possible matchgiven the value of the deviation. A maximum value of deviation of 0.92was used. This value was determined from the Matsuoka-Nakai materialmodel (a version of the MBP model). As shown in FIG. 28, a deviation of0.5 gave the closest match, hence this value was used in all subsequentsimulations.

Variation of the exponent (m) showed that the exponent value onlyinfluenced the amount of sand produced, but it did not change the pointat which massive sanding occurred. Larger values of the exponent reducedthe overall volume of sand produced (FIG. 29). Varying the value of thecritical plastic strain showed that this constant governed the onset ofmassive sanding, where larger values of the exponent resulted in sandingat later times (i.e., at larger confining stresses) (see FIG. 30).Reducing the value of the critical plastic strain resulted in sandingoccurring at lesser confining stresses.

FIG. 31 shows the evolution of the shape of the material failure fromsanding after the hollow cylinder test was performed. Each panel in FIG.31 is a horizontal slice taken at different depths in the core sample.From FIG. 31, it was observed that sanding was most massive at thecenter core slice of the sample (i.e., the cavity is deepest in thecenter of the core) and less towards the ends of the core sample. Thisshape was also observed in the numerical simulation runs. FIG. 32 showsthe simulated evolution of the shape of the material failure over timecomputed for the numerical simulation (each slice in FIG. 31 is takenvertically). Although initially failure was located at the rock-spacerinterface, the cavity progressed such that the deepest cavity occurredat the center of the core (see FIG. 32), i.e., the cavity became deepertowards the center of the simulated core, which matches what wasobserved in the hollow cylinder test results.

6.3.2 Cavity Growth

As can be seen from FIG. 32, at 116675 seconds, sand was still beingproduced. To determine if the sand production would stabilize, theloading at 116675 sec (i.e., a confining stress of 46 MPa and a pressuredifferential of 0.1643398 MPa) was held constant for an additional116675 seconds. However, the simulation was set to terminate if morethan 50% of the sample was produced a sand. FIG. 33A shows the sandproduction (cc) vs. time (s), and FIG. 33B shows the cavity size for thenumerical simulation, where the confining stress is constant at 46 MPaand the pressure differential is constant at 0.1643398 MPa from 116675seconds onwards. A continuous 0.01 mm sized grid was used. In FIG. B thegrey area indicates the location of the cavity and dark area indicatesintact rock (or spacer). At 126618.836995 seconds, 50% of the sample wasproduced as sand and the simulation terminated (see FIGS. 33A and 33B).In the volume calculation, the volume of the spacers was counted also.As shown in FIG. 33B, the cavity was very deep in the center andpenetrated almost the complete simulated rock sample. From this, it wasconcluded that a confining stress of 46 MPa and a pressure differentialof 0.1643398 MPa resulted in unstable sanding.

To determine if a stable cavity could be developed and sustained, theloading was changed. The material constants used were identical to thosein FIG. 25, i.e., deviation=0.5, beta=0, critical plastic strainvalue=0.017 and the exponent (m)=2. A similar loading procedure asoutlined in connection with FIG. 32 was used. Specifically, the loadingreached at 113230 sec, i.e., a confining stress of 44 MPa and a pressuredifferential of 1.255317 MPa, was held constant for an additional 120120seconds (max simulation time 233350 s). This simulation showed that amaximum of 10.2966 cc of sand was produced, and that the cavitystabilized (FIG. 34). This observation indicated that conditions existedunder which sand production could be expected to some extent, but whichshould not lead to catastrophic failure (due to massive sanding).

6.4 Summary

The results showed that both the onset of sanding and the massivesanding of the hollow cylinder sanding test could be fit uniquely usinga geomechanical model including a production function. The parametersused to fit the data are deviation, beta, critical strain, and theexponent. For the Landana brine test, it was found that the followingparameter values resulted in the best match:

Deviation=0.5

Beta=0

Critical strain value=0.017

Exponent=2

The results also indicated that conditions existed under which someamount of sand production could be expected, but which should not leadto catastrophic failure (due to massive sanding). The workflow for thefitting procedure in Section 6.3.1 was as follows:

-   -   1) Based on the total sand volume and the precision with which        the sand volume was measured, an appropriate grid size was        chosen. For this study, where a total amount of approximately 6        cc of sand was produced, a cell size of 0.01 mm was used to        obtain the desired precision.    -   2) Multiple models with varying deviation values were created        using the material constants obtained from fits to the data from        the triaxial test. In the computations to fit the test data, the        geomechanical model did not include the sand production model.    -   3) Using the results from step 2, the yield value at massive        sanding was recorded. This yield value was set equal to the        critical plastic strain limit, and the value of the exponent was        kept equal to 1. In this step, a number of runs could be made        with a range of values of critical plastic strain around the        obtained value.    -   4) Using the results from step 3, the simulation which fit the        onset of massive sanding was used to vary the value of the        exponent, where an increase in the exponent value was observed        to reduce the amount of sand produced. The simulation which most        closely fit the data from the tests was chosen to be        representative.

7. REFERENCES CITED

All references cited herein are incorporated herein by reference intheir entirety and for all purposes to the same extent as if eachindividual publication or patent or patent application was specificallyand individually indicated to be incorporated by reference in itsentirety herein for all purposes. Discussion or citation of a referenceherein will not be construed as an admission that such reference isprior art to the present invention.

8. MODIFICATIONS

Many modifications and variations of this invention can be made withoutdeparting from its spirit and scope, as will be apparent to thoseskilled in the art. The specific examples described herein are offeredby way of example only, and the invention is to be limited only by theterms of the appended claims, along with the full scope of equivalentsto which such claims are entitled.

As an illustration of the wide scope of the systems and methodsdescribed herein, the systems and methods described herein may beimplemented on many different types of processing devices by programcode comprising program instructions that are executable by the deviceprocessing subsystem. The software program instructions may includesource code, object code, machine code, or any other stored data that isoperable to cause a processing system to perform the methods andoperations described herein. Other implementations may also be used,however, such as firmware or even appropriately designed hardwareconfigured to carry out the methods and systems described herein.

The systems' and methods' data (e.g., associations, mappings, datainput, data output, intermediate data results, final data results, etc.)may be stored and implemented in one or more different types ofcomputer-implemented data stores, such as different types of storagedevices and programming constructs (e.g., RAM, ROM, Flash memory, flatfiles, databases, programming data structures, programming variables,IF-THEN (or similar type) statement constructs, etc.). It is noted thatdata structures describe formats for use in organizing and storing datain databases, programs, memory, or other computer-readable media for useby a computer program.

The systems and methods may be provided on many different types ofcomputer-readable media including computer storage mechanisms (e.g.,CD-ROM, diskette, RAM, flash memory, computer's hard drive, etc.) thatcontain instructions (e.g., software) for use in execution by aprocessor to perform the methods' operations and implement the systemsdescribed herein.

The computer components, software modules, functions, data stores anddata structures described herein may be connected directly or indirectlyto each other in order to allow the flow of data needed for theiroperations. It is also noted that a module or processor includes but isnot limited to a unit of code that performs a software operation, andcan be implemented for example as a subroutine unit of code, or as asoftware function unit of code, or as an object (as in anobject-oriented paradigm), or as an applet, or in a computer scriptlanguage, or as another type of computer code. The software componentsand/or functionality may be located on a single computer or distributedacross multiple computers depending upon the situation at hand.

1-20. (canceled)
 21. A method for use in predicting sand production froma geomechanical reservoir system, comprising: receiving data indicativeof plastic deformation and sanding of material within the geomechanicalreservoir system; computing, on a computer system, a value of one ormore hardening parameters based on a first fit of a hardening model tothe received plastic deformation data; wherein said hardening modelmodels plastic deformation behavior of said material within thegeomechanical reservoir system; computing, on a computer system, a valueof a critical plastic strain based on a second fit of the geomechanicalmodel comprising said hardening model to the received sanding data;computing, on a computer system, a value of at least one parameter of aproduction function based on a third fit of the geomechanical modelcomprising said production function to the received sanding data andusing a value of at least one value of said hardening parameters andsaid value of said critical plastic strain; wherein said productionfunction predicts sand production from said geomechanical reservoirsystem; and outputting to a user interface device, a monitor, a computerreadable storage medium, a local computer, or a computer that is part ofa network, at least one value of said one or more hardening parameters,said value of said critical plastic strain, or said value of at leastone parameter of said production function.
 22. The method of claim 21,wherein said production function is a function ƒ(x), wherein ƒ(x)=0 whenx=0, and ƒ(x)=1 when x=1; and wherein x is a function of the criticalplastic strain.
 23. The method of claim 21, wherein said value of atleast one parameter of said production function is a value of anexponent of said production function.
 24. The method of claim 23,wherein said sand production function is given by the expression:${f(x)} = \left( \frac{ɛ^{p}}{ɛ_{\lim}^{p}} \right)^{m}$ whereinx=ε^(p)/ε_(lim) ^(p), wherein ε^(p) is a plastic strain invariant,wherein ε_(lim) ^(p) is said critical plastic strain; and wherein m issaid value of said exponent.
 25. The method of claim 21, wherein saidhardening model is a modified Bigoni-Piccolroaz model.
 26. The method ofclaim 21, wherein said plastic deformation data is obtained from atriaxial test.
 27. The method of claim 21, wherein said sanding data isobtained from a hollow cylinder test.
 28. The method of claim 21,wherein said first fit of said hardening model to the received plasticdeformation data is obtained by a regression.
 29. The method of claim21, wherein: said second fit of said geomechanical model comprising saidhardening model to the received sanding data is obtained by solving asystem of partial differential equations that model the geomechanicalreservoir system; the system of partial differential equations comprisesa reservoir flow model and said geomechanical model comprising saidhardening model; wherein the system of partial differential equations iscoupled through a fully-expanded Jacobian; and the solving of the systemof partial differential equations includes solving simultaneously in asingle time step the fully-expanded Jacobian based upon the receivedsanding data.
 30. The method of claim 21, wherein: said third fit ofsaid geomechanical model comprising said production function to thereceived sanding data is obtained by solving a system of partialdifferential equations that model the geomechanical reservoir system;the system of partial differential equations comprises a reservoir flowmodel and said geomechanical model comprising said hardening model;wherein the system of partial differential equations is coupled througha fully-expanded Jacobian; and the solving of the system of partialdifferential equations includes solving simultaneously in a single timestep the fully-expanded Jacobian based upon the received sanding data.31. The method of claim 21, wherein the outputted at least one value ofsaid one or more hardening parameters, said value of said criticalplastic strain, or said value of at least one parameter of saidproduction function is visually displayed on the user interface deviceor the monitor.
 32. A method for sand production from a geomechanicalreservoir system, comprising operating said geomechanical reservoirsystem in accordance with a result of the method of claim
 21. 33. Amethod for use in predicting sand production from a geomechanicalreservoir system, comprising: receiving data indicative of physical orchemical properties associated with material within the geomechanicalreservoir system; generating sand production predictions for sanding ona computer system by solving a system of partial differential equationsthat model the geomechanical reservoir system, the system of partialdifferential equations comprises a reservoir flow model and ageomechanical model of the geomechanical reservoir system that arecoupled through a fully-expanded Jacobian, a sand production criterionis applied to said geomechanical model and the geomechanical modelcomprises a hardening model, the solving of the system of partialdifferential equations includes solving in a single time step thefully-expanded Jacobian based upon the received data; and outputting thegenerated sand production predictions to a user interface device, amonitor, a computer readable storage medium, a local computer, or acomputer that is part of a network.
 34. The method of claim 33, whereinthe solving in a single time step the fully-expanded Jacobian based uponthe received data is performed simultaneously.
 35. The method of claim33, wherein the sand production criterion is a critical value of a totalstrain, a critical value of a plastic strain invariant, or a maximumeffective stress.
 36. The method of claim 33, wherein the reservoirmodel and the geomechanical model are computed on a grid comprising aplurality of grid cells.
 37. The method of claim 33, wherein theoutputted generated sand production predictions are visually displayedon the user interface device or the monitor.
 38. A method for operatinga geomechanical reservoir system to control sand production from thegeomechanical reservoir system, comprising: computing, on a computersystem, a value of at least one operating parameter based on a fit of ageomechanical model comprising a production function that predicts sandproduction from the geomechanical reservoir system to data indicative ofphysical or chemical properties associated with materials within thereservoir system, the at least one value of the operating parameterindicating a condition for stabilized sand production from thegeomechanical reservoir system; and operating said geomechanicalreservoir system in accordance with the value of the at least oneoperating parameter.
 39. The method of claim 38, wherein the dataindicative of physical or chemical properties includes plasticdeformation and sanding of material within the geomechanical reservoirsystem.
 40. A non-transitory processor readable medium containingcomputer readable software instructions used for operating ageomechanical reservoir system to control sand production from thegeomechanical reservoir system, the software instructions comprising:receiving data indicative of plastic deformation and sanding of materialwithin the geomechanical reservoir system; computing a value of one ormore hardening parameters based on a first fit of a hardening model tothe received plastic deformation data; wherein said hardening modelmodels plastic deformation behavior of said material within thegeomechanical reservoir system; computing a value of a critical plasticstrain based on a second fit of the geomechanical model comprising saidhardening model to the received sanding data; and computing a value ofat least one parameter of a production function based on a third fit ofthe geomechanical model comprising said production function to thereceived sanding data and using a value of at least one value of saidhardening parameters and said value of said critical plastic strain;wherein said production function predicts sand production from saidgeomechanical reservoir system.